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FAST HIERARCHICAL 

TOMOGRAPHY METHODS AND APPARATUS 

This is a continuation-in-part of Provisional Application Serial No. 60/501,350, filed 
September 9, 2003, incorporated by reference in its entirety. 

This invention relates to tomography, and more particularly, to methods and apparatus 
for creating pixel images from projections. 

BACKGROUND OF THE INVENTION 

Tomographic reconstruction is a well-known technique underlying nearly all of the di- 
agnostic imaging modalities including x-ray computed tomography (CT), positron emission 
tomography (PET), singly photon emission tomography (SPECT), and certain acquisition 
methods for magnetic resonance imaging (MRI). It also'finds application in manufactur- 
ing for nondestructive evaluation (NDE), for security scanning, in synthetic aperture radar 
(SAR), radio astronomy, geophysics and other areas. 

There are several main formats of tomographic data: (i) parallel-beam, in which the 
line-integrals are performed along sets of parallel lines; (ii) divergent-beam, in which the 
line-integrals are performed along sets of lines that diverge as a fan or a cone; and(iii) curved, 
in which the integrals are performed along sets of curves, such as circles, ellipses, or other 
closed or open curves. One problem of tomographic reconstruction is to reconstruct a 2D 
or 3D image from a set of its line-integral projections. Another problem of tomographic 
reconstruction is to reconstruct a 3D image from a set of its surface-integral projections, that 
is, its integrals on a family of surfaces. For example, the 3D Radon transform involves in- 
tegrals of the image on a family of 2D planes of various orientations and distances from the 
origin. Some of the problems of tomographic reconstruction, and some of the reconstruc- 
tion methods, are described in standard references such as F. Natterer, The Mathematics of 
Computerized Tomography. Chichester: John Wiley, 1986; F. Natterer and F. Wubbeling, 
Mathematical Methods in Image Reconstruction. Philadelphia: Society for Industrial and 
Applied Mathematics, 2001; A.C. Kak and M. Soaney, Principles of Computerized Tomo- 
graphic Imaging. New York: IEEE Press, 1988; and S.R. Deans, The Radon Transform and 
Some of Its Applications. New York: Wiley, 1983. 
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The method of choice for tomographic reconstruction is filtered backprojection (FBP) 
or convolution backprojection (CBP), which use an unweighted (in the parallel-beam or 
Radon Transform cases)or a weighted (in most other cases) backprojection step. This step is 
the computational bottleneck in the technique, with computational requirements scaling as 
iV 3 for an N x TV-pixel image in 2D, and at least as N 4 for an N x N x iV-voxel image in 3D. 
Thus, doubling the image resolution from N to 2N results in roughly an 8-fold (or 16-fold, 
in 3D) increase in computation. While computers have become much faster, with the advent 
of new technologies capable of collecting ever larger quantities of data in real time (e.g., 
cardiac imaging with multi-row detectors, interventional imaging), and the proliferation of 
3D acquisition geometries, there is a growing need for fast reconstruction techniques. Fast 
reconstruction can either speed up the image formation process, reduce the cost of a special- 
purpose image reconstruction computer, or both. 

The dual operation of backprojection is reprojection, which is the process of comput- 
ing the projections of an electronically stored image. This process, too, plays a fundamen- 
tal role in tomographic reconstruction. A combination of backprojection and reprojection 
can also be used to construct fast reconstruction algorithms for the long object problem in 
the helical cone-beam geometry, which is key to practical 3D imaging of human subjects. 
Furthermore, in various applications it is advantageous or even necessary to use iterative re- 
construction algorithms, in which both backprojection and reprojection steps are performed 
several times for the reconstruction of a single image. Speeding up the backprojection and 
reprojection steps will determine the economic feasibility of such iterative methods. 

Several methods have been proposed over the years to speed up reconstruction. For 
example, Brandt et al. U.S. Patent No. 5,778,038 describes a method for 2D parallel-beam 
tomography using a multilevel decomposition, producing at each stage an image covering 
the entire field-of-view, with increasing resolution. Nillson et al. U.S. Patent No. 6,151,377 
disclose other hierarchical backprojection methods. While these systems may have merit, 
there is still a need for methods and apparatus that produce more accurate images, and offer 
more flexibility between accuracy and speed. 

Accordingly, one object of this invention is to provide new and improved methods and 
apparatus for computed tomography (CT) scanning. 

Another object is to provide new and improved methods and apparatus for CT scanning 
that produce more accurate images, and offer more flexibility between accuracy and speed. 



2 



WO 2005/024722 



PCTYUS2004/029857 



SUMMARY OF THE INVENTION 

These objects are achieved or exceeded by the present invention. Pixel images are cre- 
ated from projections by backprojecting selected projections to produce intermediate images, 
and performing digital image coordinate transformations and/or resampling on intermediate 
images. The digital image coordinate transformations are chosen to account for view angles 
of the constituent projections of the intermediate images and for their Fourier characteris- 
tics, so that the intermediate images may be accurately represented by sparse samples. The 
resulting intermediate images are aggregated into subsets, and this process is repeated in a 
recursive manner until sufficient projections and intermediate images have been processed 
and aggregated to form the pixel image. 

Digital image coordinate transformation can include rotation, shearing, stretching, 
contractions, etc. Reampling can include up-sampling, down-sampling and the like. 

Projections can be created from a pixel image by performing digital image coordinate 
transformation and/or resampling and/or decimation and re-projecting the final intermediate 
image. 

BRIEF DESCRIPTION OF THE DRAWINGS 

The above mentioned and other features of this invention and the manner of obtaining 
them will become more apparent, and the invention itself will be best understood by reference 
to the following description of an embodiment of the invention taken in conjunction with the 
accompanying drawings, in which: 

Fig. 1 is a block diagram of apparatus used for the present invention; 

Figs. 2A, 2B and 2C are diagrams of sampling patterns used in some embodiments of 
the present invention; 

Figs. 3A, 3B and 3C are additional sampling patterns used in some embodiments of 
the present invention; 

Fig. 4 is a diagram illustrating a known method of backprojection; 

Fig. 5A is a diagram illustrating an algorithm for one embodiment of the present 
invention; 
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Fig. 5B is a diagram illustrating the manner in which intermediate images are gener- 
ated in the embodiment of Fig. 5A; 

Fig. 6 is a diagram illustrating Fourier characteristics used to produce the intermediate 
images of Fig. 5 A; 

Figs. 7 A, 7B and 7C are diagrams showing Fourier supports of intermediate images 
for the backprojection algorithm illustrated in Fig. 5A, when the coordinate transformation 
is a digital image rotation; 

Fig. 8 is a diagram illustrating an algorithm used in another embodiment of the present 
invention; 

Fig. 9 is a diagram showing the evolution of the spectral support in the algorithm of 
Fig. 8, the blocks (1...9) corresponding to the corresponding points in Fig. 8; 

Fig. 1 OA is a diagram describing an algorithm for embodiment of the present inven- 
tion; 

Fig. 10B illustrates the image coordinate transformation used in the embodiment of 
Fig. 10A; 

Fig. 11A is a diagram illustrating shear scale backprojection, and Fig. 1 IB is a dia- 
gram illustrating hierarchical shear scale backprojection; 

Figs. 12A and 12B are diagrams showing the effect of image shearing on the spectral 
support of intermediate images; 

Fig. 13 is a diagram illustrating an algorithm for another embodiment of the present 
invention; 

Fig. 14 is an algorithm for finding optimal shear factors; 

Fig. 15 illustrates an algorithm for another embodiment of the present invention; 

Fig. 16 illustrates an algorithm for still another embodiment of the present invention; 

Fig. 17 is a diagram which illustrates common fan beam geometry with a circular 
scanning trajectory; 

Fig. 18 illustrates an algorithm for another embodiment of the present invention; 
Fig. 19 illustrates weighting of intermediate images in the algorithm described in Fig. 
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18; 

Figs. 20A, 20B and 20C illustrate sampling points for the second hierarchical level of 
Fig. 18; 

Figs. 21A and 21B are diagrams of sampling patterns used in the algorithm of Fig. 18; 

Fig. 22 is a diagram showing original intersection points obtained using the method 
illustrated in Figs. 20A-20C; 

Figs. 23A, 23B, 23C and 23D illustrate sampling points for rotation for up-sampling 
used in Fig. 18; 

Fig. 24 illustrates local spectral support at a point of an intermediate image generated 
by the algorithm of Fig. 18; 

Fig. 25 is a diagram of nonuniform sampling patterns used in the algorithm of Fig. 18; 

Figs. 26 A and 26B illustrate sampling patterns for resampling and coordinate trans- 
formation; 

Figs. 27A and 27B illustrate alternative sampling schemes which can be used in the 
present invention; 

Fig. 28 includes two diagrams of divergent beams used in the present invention; 

Fig. 29A is a diagram showing a conebeam, and Fig. 29B is a diagram illustrating 
resampling; 

Fig. 30 is a diagram of an algorithm used for resampling projections; 

Fig. 31 is another algorithm used for resampling in the embodiment of the present 
invention; 

Fig. 32 is a diagram of diagram of an algorithm used for fast hierarchical reprojection; 

Fig. 33 is a diagram of another algorithm for fast hierarchical reprojection; 

Fig. 34 is a graph showing the results of experiments using the present invention; 

Fig. 35 includes sample images generated with the present invention; 

Fig. 36A is a display of a reconstructed image using the conventional algorithm, and 
Fig. 36B shows a result obtained with the fast algorithms of the present invention; and 
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Figs. 37 A and 37B are diagrams of point spread functions of algorithms comparable 
to conventional algorithms, and Fig. 37C displays a point spread function of a fast algorithm 
of the present invention. 



DETAILED DESCRIPTION 



Symbols and Fonts 

The following system of mathematical symbols and fonts will be used to improve 
clarity. 

Functions in the space domain are denoted by small letters (e.g. /(#)), while their 
Fourier transforms are denoted by capital letters (F(uj)). 

The indices of two-variable functions are denoted variously, depending on conve- 
nience. The following three notations of function f are equivalent: f(x u x 2 ), and 

/aai). 

Continuous-domain and discrete-domain functions respectively are distinguished by 
the style of parentheses used with their indices: x 2 ) is a function of two continuous 

variables (i.e., fel 2 ), and f[m u m 2 ] is the sampled version of f{x) and is therefore a 2-D 
array. 

A linear operator and its corresponding matrix are distinguished by the font style. 
Suppose (^4 € R 2x2 ) is a matrix, then A is its associated linear operator. For example, if A 
is a coordinate transformation, g{x) = (Af) (x) = f(Ax). 

The same operator is sometimes denoted differently inside and outside block diagrams. 
While outside it may be denoted as A(a), within the block diagram it is denoted as A a . 



Overview of Hardware 
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The present invention has application in a variety of imaging apparatus, including CT 
scanners. Typical imaging apparatus 10 (Fig. 1 1) includes a scanner 12 which acquires data 
from an object such as a head, and sends raw data corresponding to line-integral projections, 
e.g., with a divergent beam geometry, 14 to a projection pre-processor 16. The projection 
pre-processor 16 applies various conversions, normalizations, and corrections to the data, as 
well as weighting and filtering, which may be shift varying. The output of the projection 
pre-processor 16 is a collection of pre-processed projections, hereinafter simply referred to 
as projections, also called sinogram! 18, which is fed to a sinogram update processor 20. The 
sinogram update processor 20 possibly modifies the input sinogram! 18, using information 
from sinogram 2 34, for example correcting for various artifacts including beam-hardening, 
or as part of a multi-step or iterative reconstruction procedure. 

The output of the sinogram update processor 20 is a sinogram 3 22, which is input 
to a fast backprojection processor 24. The fast backprojection processor 24 is generally 
a computer or special purpose hardware, or any combination thereof, of any suitable type 
programmed and/or wired to perform the algorithms described herein. 

The output of the fast backprojection processor 24 is an electronic imagei 26, which 
is input to an image conditioning processor 28. The image conditioning processor 28 per- 
forms necessary postprocessing of the electronic image, possibly including identification 
and extraction of artifact images, or images for further processing in a multi-step or iterative 
reconstruction process. 

If desired, the image conditioning processor 28 can produce an electronic image2 30 
that is fed to a fast reprojection processor 32. The fast reprojection processor 32 is generally 
a computer or special purpose hardware, or any combination thereof, of any suitable type 
programmed and/or wired to perform the algorithms described herein. If desired, this pro- 
cessor can share the use of the same computer and hardware employed by the backprojection 
processor 24. 

The output of the fast reprojection processor 32 is a sinogram 2 34, which is fed back 
into the sinogram update processor 20. The backprojection/reprojection process can continue 
until suitable results are obtained. While reprojection is not always needed, it is helpful in 
many situations. 

When the electronic imagei 26 is suitable, the image conditioning processor 28 pro- 
duces an electronic image 3 36, which is fed to a storage/analysis/display device 38. It is 
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contemplated that the electronic image 3 36 can be stored in a computer memory, and/or 
analyzed electronically for anomalies or dangerous materials, for example, and/or displayed, 
and/or printed in some viewable form. 



Overview of Backprojection and Reprojection Methods of the Present Invention 

The backprojection methods of the present invention use various techniques to create 
an image made of pixels (picture elements) and/or voxels (3D picture elements), hereinafter 
referred to collectively as pixels, which will now be introduced in a general way. 

This explanation uses terminology and processes commonly used in multi-dimensional 
signal processing, for example as described in D. Dudgeon and R. Mersereau, Multidimen- 
sional Digital Signal Processing. Englewood Cliffs: Prentice-Hall, 1983. Some terms in this 
description of the present invention are used in the following contexts. The term Sampling 
pattern refers to a set of points in space with positions defined relative to a system of coor- 
dinates. Examples of sampling patterns are seen in Figs. 2A-2C and 3A-3C. A Cartesian 
sampling pattern refers to a set of points formed by the intersection of two mutually per- 
pendicular sets of parallel lines. The term continuous image refers to a function defined on 
a coordinate system, for example, f{x, y), and f{x, y, z) are respectively 2D and 3D func- 
tions. A digital image is an array of values of a continuous image on a sampling pattern. 
More broadly, a continuous image can be represented by an array of numbers that serve 
as the coefficients in a series expansion with respect to some basis set, such as splines, of 
which the Cartesian product of zero-th order splines yields the familiar square pixel form for 
displaying digital images as continuous images. Hereinafter, this array of numbers will be 
also referred to as a digital image. All images stored in a digital computing device must be 
digital. For brevity, both digital and continuous images will be often referred to hereinafter 
simply as images, with the meaning inferred from the context. With this terminology, a pixel 
image is a digital image corresponding to a sampling pattern that is a lattice, i.e., a uniformly 
spaced, periodic pattern, usually but not necessarily Cartesian. 

One sampling pattern will be said to be sparser than another, if it yields a smaller total 
number of samples of a continuous image. Typically a sparser sampling pattern will have 
a lower sampling density. Oversampling refers to using more samples than necessary to 
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represent a continuous image at a desired accuracy. The corresponding digital image will be 
said to be oversampled. Given a digital image corresponding to a continuous image for one 
sampling pattern, the process of producing a new digital image corresponding, with a desired 
accuracy, to the same continuous image on a different sampling pattern, will be called digital 
image resampling. Up-sampling and down-sampling are special cases of resampling on a 
denser or sparser sampling pattern, respectively. Further, up-sampling or down-sampling by 
a factor of 1 involves no change in the digital image, is considered a form of up-sampling 
or down-sampling. Digital image addition refers to to a point-by-point addition of digital 
images defined with respect to the same sampling pattern, or the same basis, in the case of 
an expansion with respect to a basis. Lower dimensional digital filtering refers to digital 
filtering of a multidimensional array along only a subset of its dimensions, for example, 
separate ID filtering of each column of a 2D rectangular digital image. 

Coordinate transformation of a continuous image refers to operations such as rotation, 
shearing, stretching, or contraction. To define a digital coordinate transformation, consider 
two continuous images related by a coordinate transformation, and the corresponding digi- 
tal images representing the continuous images with respect to a common sampling pattern. 
The process of producing one digital image from the other is called digital image coordinate 
transformation. This can be accomplished by digital filtering, i.e., by discrete index opera- 
tions on the digital image. Specific examples include (but are not limited to) digital image 
rotation, digital image shearing, digital image stretching, digital image contraction, and the 
combinations of such operations. Methods for performing digital image coordinate transfor- 
mation are known, for example, as described in M. Unser, P. Thevenaz, and L. Yaroslavsky, 
Convolution-based interpolation for fast, high-quality rotation of images. IEEE Transactions 
Image Processing, Vol. 4, pp. 1371-1381, 1995. 

Some digital image coordinate transformations are illustrated in Figs. 2A-2C and 3A- 
3C. Fig. 2A shows the outline of a continuous image (a rectangle) and the sampling pattern 
for a digital image representing it. Values of the continuous image on the heavy dots are 
included in the digital image. Figs. 2-B and 2-C show the rotated and the sheared continuous 
image, respectively, on the same sampling pattern, with the heavy dots showing the values 
included in the digitally rotated/sheared version of the digital image in Fig. 2-A. 

Fig. 3 A also shows a continuous image and a sampling pattern defining a digital 
image. Fig. 3B shows the stretched continuous image by some constant factors in the x and 
y dimensions. The digitally stretched image is defined by values of the stretched continuous 
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image on the heavy dots. Fig. 3C shows the same continuous image as in Fig. 3A, but with 
a sampling pattern denser by certain stretch factors. The corresponding digital image in Fig. 
3C will be the same as in Fig. 3B. More generally, digital image stretching or contraction 
can be equivalent, for sampling patterns with some regularity, to digital image up-sampling 
or down-sampling. 

Note that the application of certain coordinate transformations, such as rotation by 0 
degrees, shearing by a shear parameter of zero, or stretching or contraction by a factor of 1, 
leave the digital image unchanged, and therefore may be either included or omitted in the 
description of a process, without changing the result. 

The Fourier transform of a continuous image allows one to determine, via sampling 
theory, the properties of sampling patterns such that the corresponding digital image repre- 
sents the continuous image to desired accuracy, as explained, for example, in D. Dudgeon 
and R. Mersereau, Multidimensional Digital Signal Processing. Englewood Cliffs: Prentice- 
Hall, 1983. Likewise, the discrete-time Fourier transform (DTFT) of a digital image allows 
one to determine what digital image coordinate transformations will produce a digital image 
that represents the transformed continuous image to a desired level of accuracy. The relevant 
properties of the Fourier transform of continuous images, and the DTFT of digital images, 
will be collectively referred to as Fourier characteristics. 

Weighted backprojection operations require weighting of each projection by a weight 
that depends on the position of the pixel generated by the backprojection. Different weights 
can be used for different projections, depending, for example, on the position of the source 
at which the projection was acquired, as explained in A.C. Kak and M. Slaney, Principles 
of Computerized Tomographic Imaging. New York: IEEE Press, 1988. As a special case, 
the weighting factor can be equal to 1, which corresponds to no weighting, but is a weight- 
ing factor. Unless specifically indicated, weighted and un-weighted backprojection will be 
collectively referred to as backprojection. 

With this background information, several embodiments of the invention will be de- 
scribed. 

The backprojection processor 28 (Fig. 1) is programmed to perform the algorithms 
used to practice the present invention. The algorithms will be discussed in detail, but will 
first be described in more general terms. Steps in the backprojection process are indicated in 
block diagrams, with reference numerals. 
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Fig. 4 illustrates rotation-based backprojection, as the sum of rotated intermediate im- 
ages formed by backprojecting individual projections gi-.gp at zero angle in step 50. The 
backprojections produces images which are subjected to coordinate transformation at 52, 
and aggregated at 54 to produce an image /. By itself, this structure is equivalent to con- 
ventional backprojection, and offers no reduction in operation count. However, it serves as 
a stepping stone to the introduction of some of the fast hierarchical backprojection methods 
of the present invention. 

Fig. 5A illustrates a hierarchical backprojection method for creating a pixel image 
/ from a plurality of projections ?i...gr. Each projection q m is backprojected at 100 to 
produce a plurality of intermediate images I lt i..Ji, P . This is the zeroth or preparatory level 
of a hierarchical backprojection. Digital image coordinate transformations are performed 
on selected intermediate images at 102. Subsets of the transformed intermediate images 
(pairs, here), are aggregated at 104 to produce aggregate intermediate images I 2 ,i, ...i"2,p/2- 
This is the first level of a hierarchical backprojection. The aggregate intermediate images 
of the first level serve as new intermediate images input to the next level of hierarchical 
backprojection. The process of applying digital image coordinate transformations to selected 
intermediate images at 106, and aggregating selected intermediate images at 108 to produce 
new intermediate images continues until all intermediate images have been processed and 
aggregated into the final image / at 1 16. 

If desired, the operations within and across some of the levels can be combined. For 
example, the zeroth and first level can be combined for some of the projections, and the 
corresponding intermediate images from the set J 2j i , --hjpp produced instead by a backpro- 
jection 1 12 at selected view angles of two or more (exactly two, for the embodiment in Fig. 
5A) selected projections q p , as shown in Fig. 5B. Alternatively, some of the initial interme- 
diate images can be produced by an equivalent process involving no explicit backprojection, 
such as a direct Fourier reconstruction method, as described in R Natterer and F. Wubbeling, 
Mathematical Methods and Image Reconstruction. Philadelphia: Society for Industrial and 
Applied Mathematics, 2001. 

The parameters of the various digital image coordinate transformations are chosen to 
account for the view angles of the constituent projections of the intermediate images, and for 
the Fourier characteristics of the intermediate images, so that the aggregates of the interme- 
diate images may be accurately represented by sparse samples, as will be explained. These 
Fourier characteristics focus on the essential spectral support of the intermediate images, i.e. 
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the region in the Fourier (frequency) domain, where the Fourier transform of the intermediate 
image is significantly different from zero. Sampling theory teaches that the spectral support 
of a continuous image determines the nature of the sampling patterns that produce digital im- 
ages that represent the underlying continuous image, and from which the continuous image 
can be reliably reconstructed. In particular, Fig. 6 shows the typical wedge-shaped spectral 
support of intermediate image in the hierarchical algorithm. 

Figs. 7(A) and 7(B) show Fourier supports of intermediate images in the binary hierar- 
chical backprojection algorithm illustrated in Fig. 5(A), when the coordinate transformation 
is chosen to be a digital image rotation. Fig. 7(A) shows the Fourier-domain support of 
the virtual intermediate images J. A virtual image ii, m is composed of the backprojection of 
the projections that are included in the corresponding image ij >m . Fig. 7(B) shows that by 
choosing the parameters of the coordinate transformation to account for the view angles of 
the constituent projections of the intermediate images, and for the Fourier characteristics of 
the intermediate images, the vertical bandwidth (the height of the broken-line rectangle) of 
the intermediate images can be minimized, allowing sparse sampling and reducing compu- 
tational requirements. Fig. 7(C) shows the space-domain sampling scheme of I^mJi-i^m 
and /i_ 1)2m _x. The sampling points are at the intersections of the dotted lines in the space 
domain. 

Fig. 8 describes another embodiment of the present invention, which performs a 
ternary aggregation of intermediate images, with P = 2 x 3 L projections. For concrete- 
ness of illustration, the case of P = 18 projections is shown. As in Fig. 5(A), projec- 
tions qe 1 yQd 2 7 -~Q0i8 816 backprojected at 100 to produce a plurality of intermediate images 
^1,1— ^i,i8- This is the zeroth level of a hierarchical backprojection. 

The projections are grouped in 3's, and selected projections (two of the three in Fig. 
8) in each group are subjected to digital image coordinate transformation at 102. Subsets 
(triplets, in Fig. 8) of the transformed intermediate images Jf, are aggregated at 104 to 
produce aggregate intermediate images J^u —i ^2,6- * s ^ l eve * °f hierarchical 

backprojection (I = 1). 

In the second level of hierarchical backprojection (I = 2) selected aggregate interme- 
diate images ii? m undergo a coordinate transformation composed of stretching or upsampling 
along the y coordinate at 106, and another coordinate transformation /C 2 , m on selected ag- 
gregate intermediate images at 108. Here, also, the transformed intermediate images are 
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aggregated in groups of three at 1 10 to produce intermediate images J^, I$ i2 - Selected inter- 
mediate images from level two are subjected to digital image coordinate transformations at 
1 12 and aggregated at 1 14 to produce the next level intermediate images Jf fl . This process is 
repeated to the extent necessary to produce the image / at 116, depending on the number of 
projections. The digital image coordinate transformation denoted by K at the last level 1 12 is 
a digital image rotation, and those at the preceding levels 102 and 108 can also be chosen to 
be digital image rotations. Here too, the parameters of the various digital image coordinate 
transformation are chosen to account for the view angles of the constituent projections of the 
intermediate images, and for the Fourier characteristics of the intermediate images, so that 
the aggregates of the intermediate images may be accurately represented by sparse samples, 
as will be explained. 

Fig. 9 shows the evolution of the spectral support in the ternary backprojection al- 
gorithm of Fig. 8. The numbers (l),.-,(9) correspond to the corresponding points in the 
block diagram of the algorithm in Fig. 8. The spectra shown are the discrete-time Fourier 
transforms (DTFTs) of the digital images J/ >m . 

Fig. 10(A) describes an algorithm for another embodiment of the present invention, 
using ternary two-shear-based hierarchical backprojection. The embodiment of Fig. 10(A) 
is similar to the embodiment of Fig. 8, and the reference numerals from Fig. 8 are used 
where appropriate. As in Fig. 8, P = 2 x 3 L projections, and the case of L = 2 is shown. 
The embodiment of Fig. 10(a) differs from that described in Fig. 8 in two respects. First, an 
additional upsampling step along the x coordinate is included at 101 in the first level coordi- 
nate transformations. Second, at least some of the digital image coordinate transformations, 
denoted by K at 102 and 108, are composed of a sequence of two digital image shear op- 
erations, the first (120) along the y coordinate (122), the second along the x coordinate, as 
shown in Fig. 10(b). 

Fig. 11A illustrates shear-scale backprojection, and Fig. 11B illustrates an equivalent 
hierarchical shear-scale backprojection. In Fig. 11A, the plurality of projections q u ? 4 > 
are backprojected at 140 and subjected to a shear-scale coordinate transformation at 142. 
The resulting intermediate images are aggregated at 144. 

In Fig. 11B, the projections q u g 4 are backprojected at 146, subjected to a shear- 
scale coordinate transformation at 148, and aggregated in subsets at 150. The intermediate 
images produced by the aggregation are subjected to a shear-scale transformation again at 
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152, and the resulting intermediate images are aggregated at 154. This process continues in 
a recursive manner until all of the projections and intermediate images have been processed 
and aggregated to form an image /. Here too, the parameters of the digital shear transforma- 
tions and sampling of intermediate images are selected to account for the view angles and the 
Fourier characteristics, so as to reduce the sampling requirements and required computation. 

Fig. 12 shows the effect of image shearing on the spectral support of intermediate 
images. Fig. 12A shows the spectral support of a certain subset of the projections as they 
should appear in the final image. Fig. 12(B) shows the spectral support of the intermediate 
image composed of the same projections, with coordinate transformation parameters chosen 
to minimize the highest radian frequency in the vertical direction. 

Fig. 13 illustrates an algorithm for another embodiment of the present invention, a 
ternary hierarchical shear-scale backprojection. Projections qeu q$ 1 are backprojected at 
160 and subjected to a shear-scale digital image coordinate transformation at 162. Subsets 
of the resulting images are aggregated at 164, and the resulting intermediate images are sub- 
jected to up-sampling at 166 and 168 digital shear coordinate transformations at 168. Subsets 
of those images are aggregated at 170, and selected resulting images are subjected to addi- 
tional coordinate transformation at 172. This process continues until all of the projections 
and intermediate images have been processed and aggregated at 174, to produce the image /. 
The final coordinate transformation shown here at 172, only involves rearrangement 

of pixels, when the sampling pattern for the pixel image / is Cartesian. 

Fig. 14 is an algorithm for finding parameters of coordinate transformations for previ- 
ously described embodiments, using the properties of the Fourier properties of intermediate 
images shown in Fig. 12. 

Fig. 15 illustrates an algorithm for another embodiment of the present invention, over- 
sampled ternary two-shear hierarchical backprojection. The embodiment of Fig. 15 is similar 
to the embodiment of Fig. 10A, and the reference numerals from Fig. 10A are used where 
appropriate. In addition to the steps of Fig. 10A, however, the embodiment of Fig. 15 
includes a downsampling step 109 in the one before the last, which is the second level in 
Fig. 15. Here too, the parameters of the various digital image coordinate transformation are 
chosen to account for the view angles of the constituent projections of the intermediate im- 
ages, and for the Fourier characteristics of the intermediate images, so that the aggregates of 
the intermediate images may be accurately represented by sparse samples. However, certain 
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degrees of oversampling are used to improve the accuracy of subsequent processing, as will 
be explained. If desired, for improved computational efficiency the downsampling step 109 
can be combined with the second, x-coordinate shear comprising the two-shear digital image 
coordinate transformation 108 (shown in Fig. 10B) producing a shear-scale transformation, 
so that the processes 108 and 109 are together replaced by the process shown in Fig. 15(B). 

Fig. 16 illustrates an algorithm for another embodiment of the present invention, over- 
sampled ternary hierarchical shear-scale backprojection. The embodiment of Fig. 16 is 
similar to the embodiment of Fig. 13, and reference numerals from Fig. 13 are used where 
appropriate. However, Fig. 16 includes additional steps of upsampling 161 in the intermedi- 
ate levels, and downsampling 169 in the level before the last level, in which the intermediate 
images are upsampled and downsampled, respectively. Similarly to the embodiment of Fig. 
15, in the embodiment of Fig. 16 when downsampling step 169 along the x coordinate 
ix Uijm follows an x-shear step 168 S£' ,m , the two can be combined for computational effi- 
ciency into a single digital image shear-scale. 

Non-cartesian Sampling Schemes For Fast Hierarchical Backprojection Algorithms 

The intermediate images in the different embodiments of hierarchical backprojection 
illustrated in Figures 8,10, and 13 have a peculiar spectral support amenable to efficient non- 
cartesian sampling. In particular, the underlying continuous domain image Jj jTn in the I th 
level occupies a wedge in Fourier space, as seen for example in Figures 6, 7, 9, and 12. Multi- 
dimensional sampling theory tells us that for images with a spectral support such as this, 
sampling on a cartesian grid is less efficient than sampling on an appropriate non-cartesian 
grid. Non-cartesian sampling can improve sampling efficiency by packing the copies of the 
baseband spectrum more tightly in the Fourier plane. For an explanation of 2D sampling see 
[?]. In particular, a quincunx sampling scheme reduces the sizes of intermediate images, and 
therefore the computational costs of the algorithm, by a factor of almost 2. Digital image 
coordinate transformation on periodic non-Cartesian sampling patterns can be executed ef- 
ficiently using one dimensional shift-invariant filters. Likewise, all the methods previously 
described for selection of the parameters of digital image transformations apply as well in the 
case of such sampling patterns. Therefore, all the embodiments previously describe extend 
to embodiments that use periodic non-Cartesian sampling patterns. 

Variants of the embodiments of the present invention already described are applicable 
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to 3D backprojection of a variety of forms of 3D projections, including the X-ray transform 
that arises for example in Positron Emission Tomography, and the 3D Radon transform, 
which arises in magnetic resonance imaging. 

The three-dimensional (3D) X-ray transform is a collection of integrals along sets of 
parallel lines, at various orientations, in 3D. Each 3D X-ray projection is a two-dimensional 
function that can be characterized by the 3D angle at which the lines are oriented. The block- 
diagram for hierarchical backprojection for 3D X-ray data is similar to the ones previously 
described, such as Figures 8, 10, 13, 15 or 16. The intermediate images in this case are 
three-dimensional, not two-dimensional as in the previously described examples. Each in- 
termediate image is sampled on a 3D sampling pattern that is sparsest in a direction that is 
an average of the 3D angles of the constituent projections. As the algorithm progresses the 
density of samples along this sparse-sampling (slow) direction increases to accommodate 
the increasing bandwidth in- that direction as explained by the Fourier analysis of 3D X-ray 
projections. Consequently at every stage in the algorithm, before the images are aggregated, 
each has to be upsampled along this slow direction. The extra dimension available in the 3D 
embodiment also provides more coordinate transformations available for use in the various 
steps in the algorithm, such as rotations in 3D. As in the 2D case, the parameters of these 
digital image coordinate transformations can be chosen to account for the constituent view 
angles and for the Fourier characteristics of the intermediate images. These digital image co- 
ordinate transformations can be decomposed into a sequence of one-dimensional operations, 
such as shears and shear-scales, as previously described. As in the 2D case, oversampling in 
any subset of the 3 dimensions may be enforced to improve image quality. 

A 3D radon transform projection is a one-dimensional function — a collection of in- 
tegrals along sets of parallel 2D planes, parameterized by the displacement of the plane from 
the origin. The view-angle of the projection is that of the 3D orientation angle of a vector 
perpendicular to the set of planes. The block-diagram of the hierarchical backprojection of 
3D radon projections is as in Figures 8, 10, 13, 15 or 16. In the first level the projections 
are Radon backprojected onto the 3D image domain. These images are constant along two 
dimensions, and therefore need be sampled only on the direction perpendicular to the two 
constant directions. When groups of 3D radon projections are combined, the bandwidth of 
the aggregate image may increase in one or two dimensions, depending on the view-angles of 
these constituent projections. It is therefore necessary to upsample the intermediate image, 
possibly in two dimensions, before coordinate transforming it and adding to other intermedi- 
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ate images. As in the 2D case, the coordinate transformations may be performed separably, 
may be combined with the upsampling operation, and oversampling may be enforced on the 
intermediate images. 

In the various embodiments of the present invention, digital image coordinate trans- 
formations and downsampling or upsampling operations may be performed by a sequence of 
lower (one) dimensional linear digital filtering operations. Furthermore, when the sampling 
patterns used have some periodicity, these digital filters can be shift invariant, as will be de- 
scribed in more detail. For computation efficiency, the digital filters can be implemented us- 
ing recursive filter structures, or using an FFT, as is known. One way to determine preferred 
values for the digital niters is using the theory of spline interpolation, as explained in M. 
Unser, A. Aldroubi, and M. Eden, Fast B-spline transforms for continuous image represen- 
tation and interpolation, IEEE Transactions on Pattern Analysis and Machine Intelligence, 
Vol. 13, pp. 277-285, 1991; M. Unser, A. Aldroubi, andM. Eden, B-spline signal process- 
ing: Part I-theory, IEEE Transactions Signal Processing, Vol. 41, pp. 821-832, 1993; and 
M. Unser, A. Aldroubi, and M. Eden, B-spline signal processing: Part II - efficient design 
and applications, IEEE Transactions Signal Processing, Vol. 41, pp. 834-848, 1993. 

The hierarchical backprojection methods of the present invention are applicable to 
a wide range of tomographic imaging geometries, such as divergent beam geometries, in- 
cluding fan-beam and cone-beam, with arbitrary source trajectories. The common fan-beam 
geometry with a circular scanning trajectory is depicted in Fig. 17. The ray source moves 
on a circular trajectory of radius D around an image of radius R. A fan-beam projection at 
source angle £ corresponds to line integral measurements along a set of rays parametrized 
by fan angle 7. Tsr is the distance along the ray of the source to the closest edge of the 
image-disc and T END is the distance of the source to the farthest edge of the disc. 

Fig. 18 illustrates an algorithm for an embodiment of the present invention, hierarchi- 
cal ternary rotation-based backprojection, applicable to fanbeam weighted backprojection. 
The embodiment of Fig. 18 is similar to the embodiment of Fig. 8, and the reference numer- 
als from Fig. 8 are used where appropriate. The embodiment of Fig. 18 differs from that 
described in Fig. 8 in several respects. First, P = 4xS L projections, and the case of L = 2 
is shown. Second, at the zeroth level, the initial intermediate images l£ m are produced from 
the fanbeam projections by weighted backprojection at 99, denoted here by W. Third, at 
the last stage, four rather than two intermediate images are aggregated. Fourth, the sampling 
patterns used in most of the early levels of the hierarchy are preferably chosen to be non- 



17 



WO 2005/024722 



PCT/US2004/029857 



Cartesian, as will be explained. Here too, the last digital image coordinate transformations 
only require reordering of image pixels, if a Cartesian sampling pattern is used for the in- 
termediate images l£ m in the one-before-last level. Also, as in Fig. 5(B), the operations in 
the zeroth and first level of the hierarchy can be combined, so that the intermediate images 
l£ m are produced by a direct weighted fanbeam backprojection of their three constituent 
projections, or by other means. 

In a selected number of levels, it is beneficial to modify the embodiment of Fig. 18, 
by including additional weighting steps before and after the cascade of upsampling step and 
digital image rotation (steps 106 and 108 in Figure 18), as shown in Fig. 19. The intermediate 
image 7 2 ,m is weighted by spatially varying weights at 180 and 182, respectively before and 
after steps 106 and 108. As will be explained, this weighting can be used to reduce the 
sampling requirements of the intermediate images, and thus reduce the computation. 

0 

The sampling scheme of the intermediate images afects the performance of the algo- 
rithm. The desired sampling scheme is one that uses the fewest samples without losing image 
information. Here too, the parameters of the various digital image coordinate transformations 
and resampling operations can be chosen to account for the view angles of the constituent 
projections of the intermediate images, and for the Fourier characteristics of the intermediate 
images, so that the aggregates of the intermediate images may be accurately represented by 
sparse samples, as will be explained. An alternative method for choosing these parameters 
is based on intersection of particular sets of rays or curves, as will be described. 

Figs. 20A-20C illustrate the progression of intermediate images through the levels 
of the ternary hierarchical algorithm described in Fig. 18, for the case of source angles 
uniformly spaced at intervals A^. The fans of the projections that make up intermediate 
images in the algorithms are shown. At the zeroth level, each intermediate image If is 
made up of a single projection with fan oriented at /3 = 0, shown in Fig. 20(A). After the 
first level of the recursion, each intermediate image is made up of three projections, with fans 
as shown in Fig. 20(B) After the second level, each image is made up of nine backprojected 
fans , as shown in Fig. 20(C). 

The intersection-based method for choosing coordinate transformations and sampling/ 
resampling patterns in the algorithm is illustrated in Figs. 21A amd 21B, for intermediate 
image 7 3>m , with constituent fans as shown in Fig. 20(C). The sampling pattern for J 3jm 
is made up of points that lie on the central constituent fan, which coincides with the one 
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shown in Fig. 20(A). As shown in Figs. 21(A) and 21(B), the sampling points for J 3j m are 
determined by the intersection of half of the central fan with the extremal constituent fans on 
the respective side of the central fan. 

It is advantageous to modify the resulting sampling pattern by applying two con- 
straints, to improve the accuracy of the backprojection, and reduce the computation require- 
ments. First, the density of samples along the rays is limited not to exceed the desired 
sampling density of the final image. Second, the sampling pattern is forced to contain on 
each ray at least one sample on the outside of the image disk on each side. The plots in Fig. 
22 displays the position of sampling points along a particular ray of the central fan at fan 
angle y r . The sampling points lie at constant intervals in 7', where 7' is the fan angle of one 
of the extramal fans shown in Figs. 21(a)(b) and (c). The points that fall on the continuous 
curve in Fig. 22 are original intersection points obtained using the method illustrated in Figs. 
21A and 21B. The points on the broken curve are those modified using the above two con- 
straints. ?Fig. 23(D)? shows an example of a sampling pattern obtained using this modified 
intersection method. The fan shown is the central fan of the intermediate image for which 
this sampling pattern has been produced. 

As in the case of the previously described embodiment of the present invention, it 
is advantageous to decompose the resampling and digital image rotations operations into a 
sequence of one dimensional operations. The blocks marked t U in Fig. 18 represent the 
upsampling of the intermediate image sampled on a central fan onto a finer set of sampling 
points on the same fan. This can be achieved by separate upsampling along each ray of the 
fan. For the cascade of upsampling and rotation marked by t U and 7C in Fig. 18, it is 
conveniet to do the decomposition into ID operation jointly, as illustrated in Figs. 23(A)- 
23P). In each of the four panels, the two dashed circles represent the boundary of the image 
of radius D and the source trajectory on the (larger)circle of radius R. The sampling points 
are represented by small circles. The digital intermediate image with central fan (Fan 0 ) and 
the sampling pattern shown in Fig. 23(A) is to be resampled and rotated to the angle of the 
central fan (Fan d ) shown in Fig. 23(D), with the final sampling pattern shown in Fig. 23(d). 
This is accomplished in the following two steps: 

(i) Upsample the image, separately along each ray of Fan a , to the sampling pattern shown 
in 23(b), which is defined by the intersection of Fan a with Fan d . These points will 
therefore lie on Fan d , as shown in Fig. 23(c); 



19 



WO 2005/024722 



PCT7US2004/029857 



(ii) Upsample the image separately along each ray of Fan d , to the sampling pattern shown 
in Fig. 23(d). 

These steps accomplish the combined operations of resampling and rotation, using ID 
upsampling operations. 

The Fourier-analysis techniques described for the embodiments of the present inven- 
tion in the case of intermediate images sampled on periodic sampling patterns are extended 
to devise spatially varying sampling schemes in the divergent beam case. These techniques 
are general enough to be applied to projections and back-projections on arbitrary trajectories, 
in both two and three dimensions. For the case of backprojection on non-periodic systems of 
lines or curves, the concept of local spectral support replaces that of spectral support. This 
is illustrated in Fig. 24 for an intermediate image produced by the backprojection of a single 
fan shown on the left side of Fig. 24. As will be explained, the local spectral of this contin- 
uous intermediate image at the indicated point parametrized by r and 9 is the line segment 
in the Fourier domain shown on the right in Fig. 24(b). The local spectral support at a point 
of an intermediate image composed of the backprojection of multiple fanbeam projections 
is shown in Fig. 25. On the left, the position of the point is indicated, as well as the range 
of view angles of the constituent projections for the intermediate image. The bow-tie shaped 
region on the right is the corresponding local spectral support. 

This analysis of the local spectral support is used to determine local sampling require- 
ments for the intermediate image. The resulting spatially nonuniform sampling patterns 
are indicated by the dotted arcs in Figs. 26A and 26B for two typical intermediate fan- 
backprojected images. The image boundary is indicated by the circle in broken line, and 
only few points need be taken outside this boundary. This local Fourier sampling method 
may be applied directly to find sampling schemes for arbitrary projection geometries over 
lines, curves or planes over arbitrary dimensions. 

The Fourier-based method for determining the sampling patterns for resampling and 
coordinate transformation for hierarchical fan-beam backprojection can be further extended, 
as will be explained, to sampling patterns lying on lines or curves other than central fans of 
the intermediate images. Examples of such beneficial sets are illustrated in Figs. 27(A) and 
27(B). 
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General Divergent-Beam Algorithms 

The methods described for fanbeam backprojection, extend directly to other divergent- 
beam geometries, including 3D cone-beam, which is one of the most important in modem 
diagnostic radiology as will be described. Similarly to the fan-beam geometry, in which a 
source of divergent rays moves on a trajectory (i.e., a circle) around the imaged object, in 
the general divergent-beam geometry a source of divergent rays moves on a (not necessarily 
circular) trajectory in 3D space around the imaged object. The detector surface measures 
line integrals through the imaged object. 

One embodiment of a hierarchical backprojection algorithm for a general divergent- 
beam geometry can again be described by a block diagram similar to Figure 18, but mod- 
ified in the following ways. First, at the zeroth level, the divergent-beam projections are 
zero-backprojected at 99 with the appropriate divergent-beam single-view backprojection 
W, corresponding to a nominal "zero" source position, producing initial intermediate im- 
ages. Second, because the trajectory of the source is not necessarily circular, the constituent 
divergent-beams of the intermediate images may not simply rotated as in the fan-beam ge- 
ometry, but also translated, with respect to each other. The coordinate transformations K in 
18 are selected accordingly. Third, depending on the presence of symmetries in the source 
trajectory and positions, there may or may not be "free" coordinate transformations such as 
the pixel re-arrangement which replaces a ?r/2 rotation in the fan-beam algorithm. 

As in the fan-beam case, the initial intermediate images are processed hierarchicaUy by 
the algorithm. Analogously to the fan-beam case, intermediate images that are close to each 
other in position and orientation are aggregated, in order that the aggregated intermediate 
image might be sampled sparsely. The 3D sampling patterns for the intermediate images 
are determined by studying the structure of constituent weighted backprojected divergent- 
beams that are rotated and translated with respect to each other. One sampling pattern, as 
illustrated in Figure 28, is where the samples are located along the rays of a divergent- 
beam corresponding to the central constituent beam of the intermediate image. The sample- 
spacing along each ray is chosen to ensure that all the constituent divergent-beams of that 
intermediate image are sufficiently sampled along the ray. Alternatively, a more general 
way is to use the local Fourier method to find the sampling pattern, as described previously 
for the fan-beam case. Knowing how an intermediate image is composed of its constituent 
projections, the local 3D Fourier structure of every intermediate image is determined. A 3D 
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local sampling matrix function at each point of the intermediate image is found that matches 
the sampling requirements for the local Fourier support, as described in the fan-beam case. 
This matrix function is then integrated (possibly numerically) over the image domain to 
determine the position of the samples. 

A separable method of up-sampling combined with rotation and translation onto a 
new sampling beam is achieved similarly to the fan-beam case. It reduces the 3D coordinate 
transformation and resampling operations to a sequence of ID resampling operations. As 
shown in Figure 29 (a), each divergent-beam may be regarded as the intersection of a set of 
vertical planar fan-beams that are distributed in azimuthal angle, with a set of tilted planar 
fan-beams at different elevation angles . The steps of the separable coordinate transformation 
are as follows (as shown in Figure 29 (b)) 

1 . The original divergent-beam is resampled onto the set of intersection of the rays of the 
original with the vertical planes of the new divergent-beam . These points are therefore 
located on the planes shared by the final sampling points. 

2. Steps 1 and 2 from the separable resampling in the fan-beam case are performed for 
each plane separately to resample onto the final set of points. 

With a suitably efficient sampling scheme, the fast hierarchical backprojection algo- 
rithm for divergent-beam can be expected to achieve large speedups with desirable accuracy. 
As in the fan-beam case one might use a pseudo-beam that is modified (e.g., with the location 
of the vertex moving farther away from the origin) in successive levels. 

As will be evident to those skilled in the art, the methods of the present invention 
are not limited to the examples of imaging geometries or specific embodiments described 
herein. The methods are equally applicable to general problems of backprojection with other 
geometries. 

Figure 30 illustrates resampling-based backprojection, as the sum of upsampled in- 
termediate images formed by generalized backprojection of individual projections at source 
positions fi p . It is similar to the rotation-based backprojection 4, but differs from it in two 
respects. First, in the first step the p th (possibly processed) projection is subjected to a 
weighted-backprojection 184 at the source position or orientation /3 P , rather than at zero, as 
is the case in 4. For example in the case of fanbeam projections, each projection is backpro- 
jected at the orientation of its source-angle. Second, before aggregation by addition at 188, 
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each initial intermediate image undergoes an upsampling operation at 186. This is neces- 
sary, because the sampling pattern of each of these P single-projection intermediate images 
is chosen to be an efficient and sparse pattern, so will usually be different for each projec- 
tion, and often non-Cartesian. Before the intermediate images are aggregated, they need to 
be resampled onto a common and denser grid. By itself, this structure offers no reduction 
in operation count. However, it serves as a stepping stone to the introduction of some of the 
fast hierarchical backprojection methods of the present invention. 

Figure 31 is another embodiment of the present invention, a resampling-based hier- 
archical backprojection for creating a pixel image / from a plurality of projections qi...q P . 
Figure 31 is the binary hierarchical version of Figure 30. First, as in the non-hierarchical 
case, each projection is backprojected 184 at its individual orientation, onto a sampling pat- 
tern suited to that orientation, producing an intermediate digital image. This is the zeroth 
level of the hierarchy. In the first level of the hierarchy, these intermediate digital images are 
upsampled at 186 to a denser sampling pattern common to selected pairs of images. This 
upsampling will usually be to a non-Cartesian sampling pattern, but can be performed by a 
sequence of one dimensional resampling operations, as shown in Fig. 23, or in Fig. 29(b). 
Selected resulting upsampled images are aggregated pairwise at 190, producing new inter- 
mediate images. In the second level, the new intermediate images are again upsampled at 
192, and aggregated at 194, producing new intermediate images. This process continues until 
all intermediate image and projections have been processed, producing after the last aggre- 
gation step 198 the final image /. As in the previously described embodiments, operations 
can be combined within and across a level. 

The sampling patterns at each stage in the hierarchy and the parameters of the up- 
sampling operation in the embodiment shown in Figure 31 may be chosen by any of the 
previously described methods. For example, in the case of fan-beam projections, one possi- 
ble sampling pattern for a given intermediate image would lie on the points of intersections 
of two fans: the first oriented at the central constituent source position; the second oriented at 
an extremal constituent source position. Alternatively the sampling pattern and parameters 
of the upsampling steps can be determined based on the Fourier or local Fourier proper- 
ties and the view angle of projections included in the intermediate images. In particular, for 
non-periodic sampling patterns, the local Fourier method described for the fanbeam rotation- 
based algorithm can be used to find the sampling patterns: knowledge of how the projections 
combine to form the intermediate images leads to the determination of the local spectral 
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support, which is used in turn to calculate the local sampling matrix function, which when 
integrated over the image domain produces the sampling pattern for the intermediate image. 

Fig 32 is the block-diagram for fast hierarchical reprojection. Reprojection is the pro- 
cess of computing tomographic projections from a given electronic image. The reprojection 
algorithm is found by applying the process of flow graph transpositionto any block diagram 
of a a backprojection algorithm, possibly with some change in weighing operations. In the 
process, operations are replaced by their adjoint or dual. The block diagrams for reprojection 
therefore appear similar to a reversed version of the corresponding one for backprojection. 
The reprojection process described in Fig. 32 is one such embodiment of reprojection, ob- 
tained from the backprojection algorithm described in Fig. 8. 

In the first level , a copy of the input image / is preserved in the top branch of the 
diagram as a top intermediate image, and in the bottom branch the image / is rotated at 
200 by -7r/2, producing a bottom intermediate image. In the second level, in the top half 
of the diagram, the un-rotated to intermediate image is subject to three separate digital im- 
age coordinate transformations at 202 some of which may leave the image unchanged, pro- 
ducing three different top intermediate images. A similar process is applied in the bottom 
branch, producing three bottom intermediate images. Each of the top and bottom intermedi- 
ate images (six in all) then undergoes a process of decimation (low-pass filtering followed by 
downsampling) at 204, producing new intermediate images. In the instance of the embodi- 
ment illustrated in Fig 32 there are only 2 • 3 L , with L = 2 view-angles, so the third level is 
the final one in the recursive hierarchy. In the third and final level the intermediate images 
are subject to separate coordinate transformations (206)some of which may leave the image 
unchanged, producing 18 intermediate images. The last step, which is not part of the recur- 
sion is different: each intermediate image undergoes a reprojection at zero degree at 208. 
Reprojection at zero degrees is equivalent to summing the vertical lines of pixels to produce 
a one-dimensional signal. These one-dimensional signals (518) are the output projections 
produced by the algorithm. 

The parameters of the digital image coordinate transformations in the algorithm are 
chosen by the knowledge of the Fourier characteristics of the intermediate images. These 
parameters are simply related to the parameters of the corresponding backprojection algo- 
rithm. It is easy to see that since the reprojection block-diagram is a flow transposition of a 
backprojection block-diagram, every branch of the reprojection block-diagram has a corre- 
sponding branch in the backprojection block-diagram, and the coordinate transformations in 
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the corresponding branches are mathematical adjoints of each other. In the version of this 
reprojection algorithm that corresponds to the two-shear backprojection algorithm in Fig 10, 
the coordinate transformations in the second level (202) is an x-shear followed by a y-shear, 
and the coordinate transformation in the last level (206) is an x-shear,f ollowed by a y-shear, 
and a fractional decimation in x (These three operations can be reduced to a shear-scale). 
The parameters of these shears are the negative of the corresponding parameter used in the 
backprojection algorithm. The parameter of the decimation in a; is the same as that of the 
upsampling in x in the first level of the baclq>rojection algorithm. In the shear-scale version 
of this algorithm (corresponding to the shear-scale backprojection algorithm displayed in 
Figure 13), the coordinate transformations in the second level (202) are shears in x (and the 
parameter of each shear is the negative of the corresponding parameter in Figure 13). The 
coordinate transformations in the final level (206) are shear-scales. The shear-scale used in 
the reprojection algorithm is the mathematical adjoint of the corresponding shear-scale used 
in the backprojection. The parameters of the decimation factors are also the same as the 
corresponding upsampling factors in the backprojection. 

Just like in the backprojection algorithms, oversampling of the intermediate images 
in the algorithm can be enforced by first upsampling the images at the beginning of the 
algorithm and downsampling them by the same factor at the end of the algorithm. Also, 
these operations can be combined within or across a level. 

Fig. 33 is the block-diagram of a decimation-based weighted reprojection. It it the 
flow-graph transpose of Fig. 31 It shows the reprojection of the image onto a set of projec- 
tions at 18 different source angles. Initially the given image is processed along two parallel 
paths. In each path the image is subject to three parallel resamplings (210) onto three differ- 
ent sparser sampling patterns. This resampling can be performed in a separable way using 
one-dimensional decimations (low-pass filtering followed by resampling onto a sparser grid). 
The parameters of the filter are determined by the Fourier characteristics of the intermediate 
image and the desired projections. Local Fourier analysis of the desired projections is used 
in the case when projections do not line on parallel straight lines. In the final level of the al- 
gorithm each intermediate image is again subjected to three parallel resamplings (212) onto 
sparser sampling patterns. Finally a weighted projection is performed on the image to pro- 
duce the projection p e . This involves a weighted sum of the pixels of the image to produce a 
one-dimensional projection. 
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Implementations and Experimental Results 

Preferred embodiments of the present invention were implemented in C programming 
language and tested in numerical experiments on a Sun Ultra 5 workstation with 384 MB 
RAM. The test image was the Shepp-Logan head phantom(a standard test image used in 
numerical evaluations of tomographic algorithms) By varying the parameters of the algo- 
rithms a tradeoff can be made between accuracy and computational cost. Accuracy refers to 
the quality of the reconstructed image. Though visual quality is not easily quantifiable, we 
measure the error between the reconstructed image and the original from which the radon 
transform was numerically computed. The measure of error used is the relative root-mean- 
square error (RRMSE). The RRMSE in reconstructing anNxN image /[m 2 , m x ] from the 
tomographic projections of /[m 2 , m i] is calculated as follows: 



For parallel-beam data, the test image was of size 256 x 256, the number of view an- 
gles was 486, and the Shepp-Logan filter (the ideal ramp filter multiplied by a sine function) 
was used to filter individual projections. Fig. 30 displays the RRMSE error versus the run 
times for the two algorithms at various values of the oversampling parameter 7 between 0.75 
and 1.0. The two-shear algorithm is represented by the circles and the shear-scale algorithm 
by the squares. Each algorithm is run using two types of filters — a third-order (dashed line) 
and fifth-order (solid line) spline filter called MOMS 16. For each flavor of the algorithm, 
as 7 is decreased the error of the algorithm decreases and the run time increases. The plot 
points that are not connected to any other are the non-oversampled versions of the algorithms 
represented by the connected points. In comparison the run time of the conventional algo- 
rithm, using linear interpolation, is 14s and the RRMSE error of its output image is 0.0486 
(worse than the fast algorithms displayed here). 

Some sample images from the output of the algorithms, for parallel-beam data, are 
displayed in Fig. 35. Columnwise from left to right, they are output images from the con- 
ventional backprojection, the two-shear and the shear-scale algorithms. An oversampling of 
7 = 0.82 was applied to the two fast algorithms. The lower row of images displays in detail 
a section of the corresponding images in the upper row. 
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The fast hierarchical algorithm for fan-beam geometry was successfully tested on the 
512 x 512 2D shepp-logan phantom. The acquisition geometry considered was with a source 
radius D = 1.06 x N = 544, 972 source angles, and 1025 equiangular detectors. The 
regular and oversampled version of the fast algorithm was implemented using a variety of 
interpolation methods. The resulting reconstructed images were compared to a conventional 
fan-beam algorithm that used linear interpolation. In all the experiments the projections were 
filtered with the Shepp-Logan filter. 

Figs. 36(A) and 36(B) display the reconstructed images from the conventional in Fig. 
36(A) and the fast algorithms in Fig. 36(B). The result of the fast algorithm is comparable 
to that of the conventional algorithm. The point spread functions of the fast algorithms are 
comparable to that of the conventional one. Fig. 37(B) displays the PSF of the conventional 
algorithm, Fig. 37(C) displays the PSF of the fast algorithm with linear interpolation and 
7 = 0.4 oversampling. Figure 37 (c) compares slices through the x-axis of the psfs of the 
conventional algorithm, the non-oversampled and the 7 = 0.4 oversampled fast algorithms. 
The similarity of the PSFs confirms the comparable image quality in the fast algorithms. 

The sampling scheme used was the fan intersection-based method, without the en- 
hancements suggested later. In addition to the shortcomings of this scheme mentioned pre- 
viously, for reasons of simplicity of implementation, numerous sample points not needed for 
the correct operation of the algorithm were used in the embodiment tested in the experiments. 
Despite these inefficiencies, for iV = 512 and D = 1.41JV, the ratio of (data-dependant) mul- 
tiply operations in the conventional algorithm to the fast (linear) one was 6.4 and the ratio of 
addition operations was 3.0. Note that the geometry used in this experiment, with a source 
very close to the origin (D = 2.06i?) is particularly challenging for this un-enhanced imple- 
mentation, because of the high density of samples near the source vertex. In most practical 
systems, the source is further away, reducing these effects. Furthermore using the alternative 
sampling schemes discussed earlier, will lead to much higher speedup factors. 

Detailed Description of Backprojection 

and Reprojection Algorithms Used In The Present Invention 

Backprojection is the process used to create images from processed projections. Re- 
projection is the reverse process, used to compute projections from a given image. Both 
operations are used in image reconstruction from projections, as seen in Fig. 1. Conven- 
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tional backprojection will be described first, followed by a description of the backprojection 
and reprojection methods of the present invention. The description will be given first for the 
case of parallel-beam projections of a 2D image, and then for the more general 2D and 3D 
geometries. 

The classic and preferred method used to estimate an image from its projections is the 
filtered backprojection (FBP) algorithm. The FBP involves first filtering the projections with 
a so-called ramp or modified ramp filter and then backprojecting those filtered projections. 
Additional pre-processing steps may also be applied to the projections before backprojection. 
A processed projection taken at view angle 0 P , where p = 1, ....P is the projection index, will 
be denoted alternately by q(t,p), or q p or q 9p , depending on whether the dependence on the 
variable t or 0 needs to be shown explicitly. For brevity, processed projections will be called 
projections hereinafter. 

The FBP reconstruction / from a set of P projections (at the angles specified by the 
components of the length-P vector 0), is described by: 

f(xi, x 2 ) = (Bgq) (zi, x 2 ) (2) 

where 

Definition 0.1 the backprojection operator is defined by 

p 

(B§q)(x u x 2 )=^2q(xicose p + X2sin0 p ,p) • (3) 

P =i 

The OiN 2 log P) coordinate-transformation-based hierarchical backprojection algo- 
rithms of the present invention are based on a hierarchical decomposition of backprojection 
in terms of coordinate transformations. The details of several preferred embodiments are de- 
scribed with different choices of coordinate transformation, concluding with the most general 
coordinate transformations. 

In particular, the 0{N 2 log P) rotation-based hierarchical backprojection algorithms 
of the present invention are based on decomposition of backprojection in terms of the rotation 
of images. 

The rotation-by-0 operator /C(0) which maps an image f(x) to its rotated version 
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(K(0)/)(x) is defined by 

Definition 0.2 (/C(0)/)(£) = f{K 0 $), 



where 



Definition 0.3 The matrix of rotation by angle 0 in the plane is Kb — [ SJj J£ 0 ° J 



Backprojection of a ID function q (t) at a single-angle 0 is denned as (B e q) {x) = nq{x! cos 9+ 
x 2 sin 0). In particular, the zero-angle backprojection (for 9 = 0) is 

{B 0 q) (x) = q(xi cos 0 + x 2 sin 0) = q(xi) =?([io ]x) (4) 

Definition 0.4 The smgle-0 p -backprojection intermediate image is I 9p (x) = (Be p q p )(x) 

As seen in Fig. 4, the backprojection in equation (3) may be rewritten as follows. Here 
q„ refers to the projection whose view-angle is 9 P . 



p 



/(*) = E(K(-W>ep)(*) 

One embodiment of hierarchical backprojection stems from the fact that the cumula- 
tive result of several successive rotations is still a rotation. In particular, 

/C(0!)/C(0 2 )-/C(^) = /C(0 X + 0 2 + ... + 0jv) ( 6 ) 

It follows from Equation (6) that with P = 2 L for some integer L, the block diagram in Fig. 
4 can be rearranged into a hierarchical tree structure as shown in Fig. 5. The intermediate 
image in the m th branch of the I th level is denoted as I l>m . In the initial level 

h,m = B 0 qm m = l,2,...,P (7) 

intermediate images are produced by backprojecting individual projections, and in subse- 
quent levels, i.e., for 1 = 1,2, 3, log 2 P 

Jl+l,m = /C(5, l2m _i)/|,2m-l + /C(<W) 42m ( 8 ) 
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the intermediate images undergo coordinate transformations (rotations, in this embodiment) 
and are aggregated by addition. For any set of projection angles 0 i9 i = 1, P, the interme- 
diate rotation angles Sj, m , can be chosen to guarantee that the structures in Fig. 4 and 5 are 
equivalent, so that the hierarchical backprojection algorithm depicted in Fig. 5 produces the 
desired image /. 

In fact, because there are many more intermediate rotation angles (free parameters) 
than view angles (constraints), there are many degrees of freedom in the choice of the £ £ , m . 
Also, for the digital images used in practice, digital image coordinate transformations are 
used, whose accuracy and computational cost depend on the choice of sampling patterns and 
implementation. These various degrees of freedom, or parameters, are used in the present 
invention to reduce the computational requirements of the hierarchical backprojection algo- 
rithm, as will be described. 

Definition 0.5 The set Ni >m of indices of constituent projections for an intermediate image 
I, >m is the set of indices of projections for which there is a path from the input in Fig. 5 to the 
intermediate image J/ )m . 

For example, as seen in Fig. 5, iY 3> P/ 4 = {P - 3, P - 2, P - 1, P}. It is easily shown, using 
equation (8) or upon examination of the block diagram, that 

Ni, m = {2 l -\m — l) + k:k = 1,2,..., 2 1 " 1 } (9) 

Thus N l>m lists the projections that make up intermediate image I l>m . If, instead of being 
processed by the method described in Fig. 5, these same projections indexed by the set N l>m 
were directly backprojected together at their respective view-angles, this would produce the 
virtual intermediate image Ii >m : 

Definition 0.6 /,, m = T,peN,, m j «p 

Upon examination of the block diagram in Fig. 5, it can be seen that the relative angle 
between projections in an intermediate image is preserved in the final reconstructed image 
/. In other words, for all intermediate images in the algorithm, i.e. VZ, m 

h,m = K{oti, m ) J|, m for some ai, m € (-*, A (10) 
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The intermediate rotation angles 6 l)jn of the hierarchical algorithm are completely de- 
termined by ai, m as follows. It follows easily from the definition of N l>m or Equation (9) that 
7V, + i, m = Ni, 2m -i U AT I>2m , and consequently by the Definition I. 



Equations (8), (10) and (11) imply that 



5l,2m-l — <*J+l,m - a l>2 m-l ^ d ^,2m = ^I+Lm ~ ai,2m ( 12 ) 



The intermediate rotation angles (5j, m ) and correspondingly the relative angles a J>m , 
are chosen to reduce the bandwidth of the intermediate images. Images of small bandwidth 
can be represented by few samples, and consequently reduce the computational cost of the 
whole algorithm. The bandwidth of these intermediate images is determined by understand- 
ing the algorithm in the Fourier domain and tracking the evolution of the spectral support of 
the intermediate images through the hierarchy. 



Tomographic projections in the Fourier domain 

The parameters of the digital image coordinate transformations are chosen to account 
for the view angles 0 P of the selected projections, as has been described. The parameters 
are also chosen for their effect on the Fourier characteristics of the intermediate images. 
These considerations are also used to determine which intermediate images are selected for 
aggregation together. 

Key to interpreting the backprojection algorithm in the Fourier domain is the projection- 
slice theorem that relates the one-dimensional Fourier transform (^i) of a projection V»f 
with the two-dimensional Fourier transform (F 2 ) of the image /. The projection-slice theo- 
rem 4 says that 

= (:F 2 /)(u>cos0,u>sin0) 

The backprojection algorithm of Fig. 4 can be interpreted in the Fourier domain. The 
backprojection-at-zero operator produces an image whose spectral support is Umited to the 
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horizontal frequency axis w x : 

(Boq)(x u z 2 ) = q(x x ) 4» {T*B 0 q){u>u w 2 ) = 2irQ(w 1 )5(u] 2 ) (13) 

It is well known that rotating an image in the space domain results in the rotation of 
its Fourier transform by the same angle 

</(£) = (B 0 g) (K e ti) & G@) = (r 2 B 0 q)(K e a) (14) 

It follows that the function of the backprojection operator in Equation (5) is to rotate 
the spectral component of each projection to the appropriate angle in the reconstructed image 
and add them all together. 

Assuming projections with one-sided bandwidth (in the t variable) equal to 7r, the 
typical virtual continuous intermediate image Jj, m in the hierarchical algorithm therefore has 
a spectral support of a wedge in a range of angles determined by the constituent projections, 
as shown in Fig. 6. In particular, let I l>m be an intermediate image in the block diagram 
of the algorithm shown in Fig. 5(A), with Ni, m = {6,6 + 1, ...,e}, i.e. the view-angles 
of the constituent projections of this image are {9 P : p = 6, 6 + 1, e}. Because of the 
relationship in Equation (10), the spectral support of J I)m is just a rotated version of that of 
I l>m , with rotation angle ai, m - It follows that the range of angles occupied by the wedge in 
Fig. 6is[7-0,7+£] = [0b-ai, m ,9 e -oii, m ]. The bandwidths of I t , m in the first and second 
coordinate, £l s and fi/, respectively, are indicated on the bounding rectangle for the wedge, 
which is shown by broken lines in Fig. 6. Thus, the choice of a Um provides a way to control 
the bandwidths of Jj, m . 

Sampling theory dictates how such a continuous intermediate image Ii >m may be repre- 
sented by a samples on a Cartesian sampling pattern. In particular, a continuous image with 
bandwidths of Q, s and % in the first and second coordinate may be sampled on a Cartesian 
pattern aligned with the coordinate axes, with a sample spacing of A a and A/ in the first and 
second coordinates such that A a < tt/O s and A f < tt/Q/. In order to minimize the number 
of samples required to represent the continuous image, the quantity 1/(A S A/) > n 2 Cl a Q,f, 
which is proportional to the area of the bounding rectangle, should be minimized. The rota- 
tion angle that minimizes this area, and therefore minimizes the sampling requirements for 
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intermediate image Ii >m , is 

<m = % + 0 e )/2 (15) 

Selecting Parameters for Coordinate Transformations 

Equation (15) is used in Equation (12) to determine the optimum rotation parameters 
5j m for the digital image coordinate rotations in the rotation-based hierarchical backprojec- 
tion embodiment illustrated in Fig. 5. The indices b and e are determined as corresponding 
to the smallest and largest, respectively, in the set N l>m . With this choice, the intermediate 
image J i>m = /C(a? m )Ji >m , and its bandwidths in the first coordinate (slow bandwidth) and 
second coordinate {fast bandwidth) are, respectively, 

To further minimize the bandwidth, the differences 9 e - 6 b should be minimized, which 
is achieved by selecting to aggregate at each step intermediate images produced from pro- 
jections with the least maximum angular separation between their view angles. 

In addition to rotation angles and selection of which intermediate images to aggre- 
gate, it is necessary to choose appropriate sampling patterns for the digital image coordinate 
transformations of the various intermediate images in the algorithm. These are chosen us- 
ing the sampling requirements, which are in turn determined using the spectral supports and 
bandwidths of the intermediate images. 

In particular, for initial intermediate images formed by the backprojection of a single 
projection, Ji, m = B 0 qe m as in Fig. 5, the slow bandwidth £l a (Ji, m ) = 0, and a single sample 
along the x 2 coordinate suffices. The slow bandwidth will be larger and more samples along 
the x 2 coordinate will be needed for initial intermediate images formed, as in Fig. 3-b, by 
backprojecting more than one projection. 

The sum of two virtual intermediate images (7i, m = Jj-i,2m-i + Ii-i^m) and their cor- 
responding intermediate images are illustrated, in the space and frequency domains, in Figs. 
7(A)-7(C). Figs. 7(A) and 7(B) show the Fourier-domain support of I and i". TFig. 7(C) 
shows the space-domain sampling scheme of Ii, m Ji-i,2m and Ii_i,2m-i- The sampling points 
are at the intersections of the horizontal and vertical lines in the space domain. The aggre- 
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gated intermediate image Ji >m , having a greater slow bandwidth than each of its components, 
requires a denser sampling pattern. Conversely, intermediate images at earlier levels of the 
hierarchical algorithm require a sparser sampling pattern, leading to computational savings. 



Computation Cost and Savings 

The computational cost is readily estimated. Assuming equally P projection view- 
angles in [0, tt), with an angular spacing of A* = n/P, and Ni itn as specified in (9), equation 
(16) (using b = 2 l ~ l {m - 1) + 1 and e = 2 l ~ 1 (m - 1) + 2 1 " 1 ) implies that n s (Ii, m ) = 
tt sin (tt^'- 1 - 1)/(2P)) and f2/CM = tt. The size (number of points) of the digital 
image If- in the I th level is therefore 

size^J « (N/A S )(N/A f ) = {N^ 8 (I lim )M(Nnf{Ii,m)M 

= iV 2 sin(7r(2'- 1 -l)/(2P)) 

< JV 2 (7r2 J -V(2P)) = 0{N 2 2 l /P). 



Given that the cost of rotating a digital image of S pixels is 0(S) and size(/£ m ) — 
0(N 2 2 l /P), the complexity of the algorithm is 

32 ^rO(^) = E °( iv2 > = °( N lo s p ) (17) 

Z=l *=i 

For the typical P = 0(N), this is 0(iV 2 log iV), which is much more favorable scaling 
of the cost with the size of the image, than the 0(JV 3 ) of conventional backprojection. 



Improved ternary rotation-based hierarchical backprojection 

It is easy to see that the rotation of a sampled image by certain special angles — 
0,±7r/2 and it — are computationally inexpensive operations because they involve no inter- 
polation but, at most, a mere rearranging of existing pixels. The computational efficiency 
of the algorithm of the present invention can be improved by incorporating these free opera- 
tions into it. The binary algorithm (Fig. 5) can be modified such that three, not two, images 
are added at every stage. The center image of each such triplet is rotated by 0 radians — a 
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free operation. To use the free rotation by -tt/2, it is included in the final stage: one of the 
constituent images is rotated by -7r/2 and added to the other un-rotated one. This partic- 
ular combination of binary and ternary stages results in a set of view-angles of size 2x3 L 
for some integer L. Though the algorithm can be tailored to arbitrary sets and numbers of 
projections, the description and analysis of the algorithm can be simplified by assuming that 
they number exactly 2 x Z L and are uniformly distributed as follows: 

0 = _ZL(i _ 3 -i) + A e (i - 1) where i = 1, 2, 2 • Z L and A e = tt/(2 • 3 L ) (18) 
4 

Hence, the view-angles can be divided into two sets of 3 L each, one set centered around 
0 = 0 and the other centered around 0 - tt/2. The block diagram of this ternary algorithm is 
shown in Fig. 8 for the particular case of L = 2, i.e., for a set of P = 2 x 3 2 = 18 projections. 
The digital intermediate images l£ m are the sampled versions of the underlying continuous 
intermediate images Ji, m . In the block diagram, blocks marked Bo represent the zero-angle 
backprojection operator. Blocks marked K s represent the digital image rotation operator. As 
will be explained, the sampling periods and rotation angles for the digital images and for the 
digital image rotations are chosen to nfinimize computational cost. 

Selected aggregate intermediate digital images can be stretched to produce new in- 
termediate images. The blocks labelled ty U l>m in Fig. 8 represent up-sampling along 
the second (slow) coordinate, which are necessary to accommodate the increasing band- 
width of the intermediate images as the algorithm progresses. For the configuration of 
view-angles in (18) and (22), the bandwidth of the intermediate images in the I th level are 
aCTi.m) = Trsin (A^'" 1 - l)/2) and a f (Ii, m ) = tt. The normalized spectral supports at 
key points of the algorithm, numbered (1),...,(9), are displayed in Fig. 9. As the algorithm 
progresses (as I increases), the slow bandwidth increases, and the required slow-sampling 
period, A s < ir/Sl s (Ii,m) = sin ((3 ! - x - l)/2)~\ decreases. The upsampling factors U t , m in 
the various upsampling blocks can be chosen to satisfy these requirements. 

Section 0.0.0.1 describes how the separable rotation of discrete images involving two- 
shears is incorporated into the backprojection algorithm. In order to improve the quality 
, of the backprojected image, a systematic oversampling of the intermediate images is intro- 
duced. Section 0.0.0.1 describes the algorithms modified to include this oversampling. 
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Intermediate rotation angles for coordinate transformation can be chosen as follows. 

For general L, the hierarchy consists of (L + 1) levels. In the zeroth level, each projection q g 
is backprojected-at-zero (B 0 ) to produce a corresponding image B 0 q 9 . The images are then 
grouped into threes and combined to produce a third as many images. The groupings at level 
J are defined by the following relationships. In subsequent levels, I = 1, L, 

Il+l,m = lC{5l,3m-2)Il,3m-2 + JC(Sl,3m-l)Il,3m-l + ^(Sl,3m)Il,3m (19) 

and Il+2,i = /C(<5i+i,i) h+\,i + /C(<5l + i )2 )/l+i,2 ( 20 ) 

For the configuration of view angles in (18), the optimal intermediate rotation angles are 
Si, 3m -i = 0 for I = 1, 2, L , 5^+1,1 = 0, and S L+1 , 2 = -tt/2. Upon examination of Fig. 
10, it is easy to see that the indices JV l>m of the projections that constitute image I l>m are as 
follows: 

N l>m = {3 l - x (m - 1) + 1, 3'- x (m - 1) + 2, 3»" l >m} for I = 1, 2, L + 1 (21) 

By Equation (21) therefore, b = ^{m - 1) + 1 and e = ^-^m - 1) + 3'" 1 and, using 
Equations(18) and (15) a? >m = -f (1 - 3~ L ) + M*' 1 ^ ~ 1/2) - 1/2). Similarly to the 
binary it follows that for I = 1, L, 

+A03'- 1 if i = 2 

5l,3m-i = a*+i, m - <3m-i ={o if i = 1 ( 22 ) 

-Agd 1 - 1 if i = 0 

For the final stage we make use of the free rotations by 7r/2 and 0 radians. 



0.0.0.1 Coordinate transformation can be based on shearing, as follows. The digital 
image rotation of the intermediate images can be replaced by a sequence of two digital image 
shears as shown in Fig. 10(B), where shears in the x and y coordinate are defined as follows 
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Definition 0.7 The x-shear and y-shear operators are defined by 



(*(«)/)(*) 



where the x-shear and y-shear matrices are S£ = [ \ f ] and S$ = [ * ? ]. respectively. 

This alternative embodiment is derived from the two-shear factorization of the rotation 
matrix: K e = S** 9 Sz** 9tM9 S c , where S c = [ co o s9 lf Le]- This decomposition implies 
that in the case of digital images, using perfect bandlimited interpolation and a sufficiently 
bandlimited digital image as the input to the two-shear digital image coordinate transforma- 
tion, the twice-sheared image is simply the image rotated by 0, and effectively downsampled 
in a; and upsampled in y by a factor of 1/ cos 6. To correct these changes in sampling pat- 
tern would require non-integer resampling every time a two-shear operation is performed. 
Instead, by sampling the intermediate images differently, only the cumulative effects need to 
be addressed for the final image. 

A convenient choice is to correct the resampling once at the beginning of the hier- 
archy, by upsampling the initial digital intermediate images in a; and downsampling in y. 
When the initial intermediate images are produced by single-view backprojection as in the 
embodiments shown in Fig. 10, the downsampling in y is free, because it only involves 
initial backprojection onto a sampling pattern with a different y density. It is therefore not 
shown explicitly in the block diagram in Fig. 10. Upsampling in x in the first stage avoids 
aliasing problems in later stages. For each of the initial images, the cumulative effective 

downsampUng/upsampling factor 1/ II&Li cos to 61131 ima S e is calculated across 311 
two-shear transformations on the path to the final image, and the reverse operation is applied 
at the initial image as described. It follows that each of the initial images l£ p ,p = 1, P is 
resampled with different x and y sampling densities. 

So the resulting algorithm, as shown in Fig. 10(A) involves first resampling all the 
images B 0 qe p ,P = 1, P by the appropriate factors (depending on 9 P ) and then performing 
two-shear digital image coordinate transformations as shown in Fig. 10(B). The two shears 
are performed taking into consideration the changing sampling pattern of each intermediate 



image. 
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Shear-scale algorithms 

The shear-scale-based hierarchical backprojection algorithm is based on the definition 
of backprojection as the sum of single-projection images that are scaled (stretched/contracted) 
and sheared. 

Definition 0.8 The x-shear-scale operator C (a) is defined by 

(C(v)f){x) = f(Crf) 
where the x-shear-scale matrix C ff = [*o °i ]• 

So, Equation (3) may be rewritten as 

p 

/(f) = Y^[C(cose py sme p )B Q q P )(x) (23) 

The following are some results describing the effect of accumulated shears and shear- 
scales that can be easily shown. 

n«S,(a i )=5,(t«i). < 24 > 

i=l i=l 

• ■ 

C(a')C(a") = C (oio-i', a'[ + a[<4) (25) 

Similarly to Equation (24), Equation (25) can be inductively applied to prove that a 
sequence of n shear-scales is also a shear-scale. Using this property, Equation (23) can be 
rewritten in a hierarchical algorithmic structure. Fig. 1 1(A) is the block-diagram representa- 
tion of Equation(23) (for P^=4), and Fig. 1 1(B) is its hierarchical equivalent. 

The requirement of equivalence of the two structures sets up a system of equations 
between the given projection-angles 9 P and the unknown shear-scale parameters a^ m , whose 
solution always exists. 

It is computationally beneficial (in terms of operations or storage space) to use an 
integer scale-factor rather than a non-integer one. For example, it is beneficial to use x- 
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shear-scales C{a) with an integer shear-parameter ai = 1.0, which is exactly a pure x-shear 
S x {a 2 )> It is possible to design the algorithm so that only pure x-shears are used, and x- 
shear-scalings are avoided, in levels I = 2, 3, L. 

Instead of replacing each shear-scale in Equation (23) with a sequence of shear-scales, 
using Equations (24) and (25) each shear-scale is replaced with a single shear-scale followed 
by a sequence of shears: 

C(cos0 p ,sin0 p ) = S x (a n )S x {a n -i)...S x (ai)C(cos9 p ,a) 

" . n (26) 

s.t. a + cos 9 P &i = s ^ n Op 

i=l 

with free parameters 0:2, ct n . 

The resulting hierarchical configuration is displayed in Fig. 13. A convenient choice 
is to perform the x-shear-scale (C<?) in the level I = 1 and x-shears in subsequent levels 
I = 2, 3, L. So the expression for the underlying continuous-domain image in the first 
level is as follows. 

f*^ = J^C{i^)B Q q^ (27) 

For I = 2, 3, L the intermediate image I l+X>m is composed of the x-sheared versions of 
three intermediate images ({i*,3m-i}L 0 ) in ^ P^vious level: 

/i+l,m = <S x (aj,3m-2)A3m-2 + h^m-l + Sx(<*lfrn) ^,3m (28) 

And in order to use the free rotation by 0 and —ir/2 radians, in the final level: /l+2,i = 

Xt+i f i + /C(— 7r/2)Ji + i,2. 

The parameters of the coordinate transformations in the algorithm displayed in Fig. 13 
(for the various shear-scales and shears) are underdetermined. In fact, only the (nonintegral) 
scale-factors a x in the first level are completely determined by the angle of the associated 
projection. The intermediate shear-factors ({<**, m }i=2,3,...,x,) are chosen to minimize the band- 
widths of the intermediate images, in order that they may be sampled most sparsely. Given 
any such set {a z , m }, the a 2 's (shear-factors) in level I = 1 can always be chosen to ensure 
that the backprojected image / produced by the entire hierarchical algorithm of Fig. 13 is 
equal to the backprojection expression (23). In particular, the shear-scale parameters in the 
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initial level of the hierarchy depend on the set {&i im } as follows : 

cos 0p if p < y 

cos(^-f) ifp>f 

sin 0 P - cos 0 P J2t= 2 a i,n(p,i) : P<2 
sin (0 P - § ) - cos (Op - § ) Ef= 2 ^,m(p,0 : P > f 

where /z(p, Z) = [p/3 z_1 ] describes the path in the hierarchy from the p th projection in level 
1 to the root. 

By examination of the upper half p < P/2) of the block-diagram in Fig. 13, following 
the sequence of shears backwards from the last level, it can be seen that in this top half of 
the hierarchy the underlying continuous intermediate image I/, m is related to its associated 
virtual intermediate image Ii i7n by a shear i.e., 

h t m = $x ( A ,m) hm (29) 

Similarly to the derivation of Equation (12), it can be shown that the intermediate 
shear-factors depend on the { A,m} as follows : 

Oil,3m-i = A+l,m - A,3m-i for i = 0, 1, 2. (30) 

The freedom in the choice of the shear coefficients aj, m provides a freedom in the choice 
of the /3 parameters, which can be used to minimize the bandwidth of the intermediate pro- 
jections. This will reduce their sampling requirements, and improve the computational effi- 
ciency of the algorithm. 

Let the spectral support of the intermediate image I/ j3m _i be denoted by Wj j3m _i, and 
its y bandwidth (i.e., slow bandwidth) be denoted by €ly(Wifl m -i). The optimal A, 3m _i 
minimizing fi y (Wi >3m _i) is then, for i = 0, 1, 2 

PbmM = ar § rr^{Q y (S y (P)W l9 z m ^)} (31) 

If the view angles are uniformly distributed, it is found that fii^^i « A*+i,m- An 
advantageous choice is then ^* 3m _ 1 = A*+i,m> yielding a£ 3m _ a = 0, which eliminates about 
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one third of the required shear operations. 

The set {Pi, m } is determined by tracking the spectral support of the intermediate im- 
ages in the algorithm, backwards from the final to first level of the hierarchy. Following from 
the slice projection theorem, it is clear that the spectral support of Ii, m is 

Wi >m = {(r cos 9 P , r sin 9 p ):pE N l>m , r € [-tt, tt]} (32) 

For the sake of sampling requirements, the spectral support Wj >m is characterized by its end 
points denoted as Ei tm , 

Ei, m = { (cos 9 P , sin 9 P ) : p e N l>m } (33) 

and shown in Fig. 12. The upper and lower band-edges of this set E lytn are denoted by 
E+ m and E^ m respectively. So if Ni, m = {M + l,-,e}, E+ m = (cos0 e ,sin0 e ) and 
E~ = (cos 06, sin 0 b ) . Because of the ternary combination rule, the upper, middle and lower 
third sets of these spectral-support end-points are denoted as E? +1>m ,El +l)jn , and E l+hm i.e. 
Ef+x m = E l>3m -i. And following from (30) and Fourier theory about affine transformations, 
Ei, m = S y (Pi,m)Ei, m . The optimal shear factors are then 

(■ g i+l,m)« + (ffi+l.tJy £ = o,2 (34) 

The algorithm FlNDALPHAS shown in Fig. 14 finds the optimal shear-factors a? )7n 
while traversing through the hierarchy down from level L to level 2 in the shear-scale back- 
projection algorithm (Fig. 13). The inputs to the algorithm are the sets of end-points of the 
spectral support of the images at the start of the (L + l)-th level: -E^+i.m for m = 1, 2. In 
particular, in the case of uniformly distributed angles, E L +i,i = #l+i,2 = {(cos 9 P , sin 9 P ) : 
p = l,2,...,P/2}. 



Finding Upsampling Factors 

The y-upsampling factors {Ui, m } in the algorithms and different embodiments of this 
invention determine the sampling density of the intermediate images in the slow direction. 
These upsampling factors can be chosen to meet the sampling requirements of the interme- 
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diate images, while being restricted to integer values. This restriction to integers has com- 
putational advantages, and may be expanded to that of rational numbers with low-valued 
denominators for the same reason. 

The problem of choosing these factors is slightly simpler in the case of the shear- 
scale algorithm than the two-shear case. In the former case, examining Fig. 13 it is easy 
to see that the slow sampling interval of the intermediate image I i>m is Yl v ^i Uvtfyim)* 
where ii(V , Z,m) = \m/3 l '~ l ] describes the path from 7 z>m to the final node of the algorithm. 
Sampling theory requires that slow direction sampling interval be as follows : 

L 

A'/ 71 = II Uv^V^m) < ir/fly (Wi f m) (35) 
l'=l 

where Sl y (Wi, m ) is the y-bandwidth of J i>m . The computational cost of the whole algorithm, 
given a set of upsampling factors U = {Ui )7n G Z : i = 2,3, L\ m = 1, 2, 6-3 }, 
can be shown to be M 

J(U)=± £ ( L ^ ) + c (36) 

where the constants c l>m denote the relative computational cost of digital coordinate trans- 
formation applied to intermediate digital image lf >m . These constants are determined by the 
order of digital filters used, and by the particular implementation of the coordinate trans- 
formation, as will be described in the discussion of efficient implementations of coordinate 
transformations. The best set U to minimize J(U) subject to the constraint Equation (35), 
can be solved by dynamic programming or a comprehensive search. For a given set of view 
angles, the best set of upsampling factors can be precomputed and stored in a lookup table. 

In the case of the two-shear rotation algorithm, the problem is complicated slightly by 
the fact that the use of two-shear rotation causes a defacto fractional up and down sampling 
of intermediate images. In particular, the slow-sampling periods A Z)TIl in adjacent levels of 
the algorithm are related as follows: A l s +1,lmm = A' i? ' 7l /(C/|, l »«i )m ) where 



1 if m = 2, 5, 8, 2 * 3 fc - 1 (no rotation) 

1 / cos (A* * 3'" 1 ) if m = 1, 3, 4, 6, 2 * 3* (rotation by ±A* * 3 1 " 1 ) 
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So the constraint on the slow direction sampling interval becomes as follows: 

L 

AT = (C/^(i' (i ,m)«»'^M,m)) < 7T/fi y (W Iim ) (37) 



l'=l 



where ft y (Wi, m ) is the y-bandwidth of I />m . The computational cost given by Equation (36) 
now can be minimi2ed subject to the constraint in Equation (37). 



Oversampled Versions of the Algorithms 

Oversampling can be incorporated into this algorithm to improve the accuracy of the 
reconstruction. One way to do so is uniformly over all intermediate images in the algorithm: 
whenever an interpolation is performed in the algorithm, the image being operated upon is at 
least oversampled by some predetermined constant 7. In particular, the ratio of the Nyquist 
rate to the sampling frequency (in both slow and fast directions) of every intermediate image 
that is subject to a digital coordinate transformation or resampling should preferably be less 
than 7. This oversampling is preferably incorporated while modifying the algorithms of the 
present invention as little as possible. 



Oversampling in the slow direction 

In the previously described algorithms of the present invention, the sampling fre- 
quency in the slow direction is controlled by the upsampling factors Ui, m - This proves to 
be a useful tool for maintaining the oversampling condition in the oversampled versions 
of the algorithms, and therefore results in no alterations to the structures of the algorithms 
of the present invention. The upsampling by a factor of 1/7 is achieved by simply modi- 
fying the constraint on the slow-direction sampling interval by the factor 7, i.e., requiring 
A i s ,m < 77r /ft fi (4 m ) for I = 2, 3, L. 



Oversampling in the fast direction 

In the fast direction, the computationally inexpensive integer upsampling is not used to 
control the oversampling, so the algorithm is modified to involve fractional upsampling. In 
the two-shear hierarchical backprojection algorithm, such a fractional upsampling by factors 
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{Ui >m : m = 1, 2, 2.3 L } in the slow direction is already included in level I = 1. We there- 
fore simply increase the upsampling factors Ui, m to incorporate this oversampling and then 
downsample the image after the last digital coordinate transformation has been performed, 
to return the image to the desired sampling scheme (where A/ = A s = 1.0). This modi- 
fication to the two-shear algorithm therefore involves only one additional level of fractional 
resampling. This x coordinate resampling is combined with the digital x-shear in the L th 
level for improved computational efficiency. The block diagram of this algorithm is in Fig. 
15. 

In the shear-scale hierarchical backprojection algorithm, the z-upsampling operations 
are combined with the x-shear-scales (C(a) at the beginning of the algorithm) to avoid an 
extra level of resampling.. It is easy to see that z-upsampling is simply a special case of 
x-shear-scaling, with the shear factor set to 0. The combination of these two operations is 
still an x-shear-scale: C{a u a 2 )C(l/U, 0) = C{o x /U, a 2 /U). The downsampling at the end 

of the algorithm is combined with the rc-shears in L th level (S aL ,Jm=i e). to make them 

effectively shear-scales of the digital image. It can be shown that an x-shear by a followed 
by a x-downsampling by U is effectively an x-shear-scale: C(U, 0)S x {a) = C(U, 1 -I- Ua). 

The exact values of these upsampling and downsampling factors is determined both by 
the parameter 7 and the spectral structure of the intermediate images. In the fast direction, 
the oversampling condition is that the fast direction (y coordinate) sampling interval satisfy 
A l i m < 77r/ft/(ii, m ) for I = 2, 3, L. So the additional upsampling necessary is 

_ max 2 <KL £lf(Ii,m) (38) 

Consequently, the upsampling factors U l>m , of the oversampled algorithm, are modi- 
fied from Ui >m , those of the non-oversampled two-shear algorithm, as follows: 

Ui,m = vU hm and U L +l,m = V < 39 ) 

In other words, the upsampling factors in the first level are modified to ensure that 
all the intermediate images in the algorithm are oversampled according to the parameter 
7. Consequently, after the last rotation, the images have a fast-sampling interval of l/*7'» 
and, therefore, at the L th level, they have to be downsampled by 77 to return them to a unit- 
sampling scheme. 
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It must be noted here that because the input projections are assumed to be sampled ex- 
actly at the Nyquist rate, the oversampling condition in the fast direction will not be satisfied 
for the first-level images, B 0 qo p , V — 1> ■•> P- 

The block diagram of this ternary oversampled shear-scale-based algorithm is in Fig. 

16. 

Collapsing sequences of similar operations 

Whenever there is a sequence of two operations operating on the same single coordi- 
nate, they may be combined for improved computational cost and resampling accuracy. The 
following is another example, in addition to the previously mentioned ones of combining x- 
shears andx-upsampling/downsampling, orx-shear-scales withx-upsampling/downsampling. 

The non-oversampled two-shear algorithm of Fig. 10 involves a y-shear followed by 
an x-shear of four out of six of the intermediate images in the L th stage. We incorporate the 
final downsampling of the image in the oversampled version of Fig. 15 by collapsing the 
sequence of two operations — the x-downsampling followed by the x-shear of this stage — 
into a single one. This leaves the length of the cascade of interpolations unchanged from the 
non-oversampling case. 

Hierarchies of Arbitrary Radixes/Branching factors 

All these algorithms can be easily modified for the case when the set of view-angles 
is not of the form 2 * 3 L . Though the preferred embodiment is for all the branches of the 
hierarchy to involve the aggregation of triplets of intermediate images, or use rotations by 
±7r /2 or 0, arbitrary numbers of intermediate images may be combined at each stage. 

Given an arbitrary number (say M) of intermediate images at a particular level of the 
algorithm, 3 x |M/3J of them may be combined in groups of three (where [x\ is the largest 
integer less than or equal to x). If the remaining number of intermediate images is two then 
they may be aggregated as a pair to produce an image at the next level. If there is only 
a single intermediate image remaining it may be passed on, without alteration, to the next 
level of the hierarchy. 
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The branching factor of the hierarchy (the number of branches that aggregate at a node 
of the hierarchy) may be altered, not just to accommodate arbitrary numbers of view-angles, 
but also to reduce the depth of the hierarchy and thereby improve image quality. In that case 
it may be useful to aggregate images not in pairs or triplets but in larger groups. 

The previous prescriptions for the parameters of the coordinate transformations may 
be easily extended to nodes with arbitrary branching factors. For example in the rotation- 
based algorithm where I l+1>m = ££i /C(<5 t)mi )I t>m< , the relations described in Equations 
(10) and (15) still hold, and the rotation angles therefore are prescribed by 6i, mi = aj+i.m - 
on >mi . In the case of the shear-scale based algorithm where Ii+i, m = ]Ci=i ^'*( 0! J,m 4 )^i,m 4 » 
Equation (29) still holds. Extending the notation Ef tm to refer to the i th set of the M sets of 
projections being aggregated, one obtains an equation for af >m . identical to the right-hand- 
side of Equation (34). 



0.1 Hierarchical Algorithms Based On Other Image Transformations 

The back-projection equation (3) may be written in the form 

m = ]pf>^As p x) (40) 

for any matrix A g , as long as the first row of A 0 is [cose sin*]. The freedom to choose 
arbitrary values for the remaining two entries of A e allows for flexibility in the design ,of the 
coordinate transformations used in the hierarchical algorithm. Matrix Ag can be factored as 
A g = Y[f =1 A(5i) for some parameter vectors <5j that are related to 9, but are not completely 
determined, owing to the freedom in the two bottom entries of A B . This factorization can 
be used to derive a hierarchical decomposition of the backprojection equation (40) with 
a corresponding block diagram such as Fig. 5, with the coordinate transformation steps 
denoted by K(6i, m ) representing the image coordinate transformation operators defined as 
(/C(*i,m)/)(2) = f(M\m)x). 

Clearly, the specific embodiments of this invention described herein are special cases 
of this more general choice of matrix A{ B and its factorization, and the associate coordinate 
transformation. The effect of these coordinate transformations on the Fourier spectrum of 
the intermediate images is analyzed similarly to the cases already described, because the 
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effect of an affine transformation by matrix A in the spatial domain is an affine transfor- 
mation by matrix A~ T , i.e., the inverse-transpose of A, in the frequency domain. Similar 
considerations can therefore be used to select free parameters in the transformations with the 
goal of reducing the computational requirements. Thus, the class of digital image coordinate 
transformations used in the hierarchical backprojection algorithms of the present invention 
includes many other transformations in addition to those described for specific embodiments. 
Furthermore, because matrices A{5^ m ) can be factorized into a product of triangular matri- 
ces, the coordinate transformations can be performed as a cascade of single-axis coordinate 
transformations, if desired. 

0.2 Efficient and accurate implementation of Digital Image Coordinate 
Transformations and Resampling 

The accuracy and speed of the hierarchical backprojection and reprojection algorithms 
of the present invention depend on the specifics of the implementation of the various digital 
image coordinate transformations and resampling operations. Improved accuracy requires 
in general high order filters or interpolations, which usually increases the computation cost. 
The cost of high order filtering or interpolation can be reduced by decomposing the op- 
erations into lower order operations. Additional reduction in computation and/or memory 
requirements is obtained if the filters used are shift invariant. High order finite impulse re- 
sponse filters can be implemented efficiently using low-order recursive filter, or by using the 
fast Fourier transform (FFT): 

In particular, if a separable representation basis is assumed for the continuous image, 
a digital y-shear of the digital image can be achieved by fractional delays of the individual 
vertical lines in a digital image array. Likewise, a digital x-shear can be expressed as a 
fractional delay of one row at a time. In turn, a fractional delay of a ID signal can be 
accomplished using a shift-invariant filter. Similar decompositions are known for 3D images 
and shear operations. 

Resampling a digital image can also be usually performed using lower dimensional 
operations, which can often be shift invariant. For example, image upsampling along one 
coordinate by an integer factor U may be decomposed into U different, computationally 
efficient, fractional delays. This is essentially, the so-called polyphase decomposition, well- 
known in digital signal processing. Rational, but non-integer resampling along one coordi- 
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nate can be decomposed into a cascade of integer down and upsampling, each of which is 
efficiently performed. 

More general digital image resampling can also be decomposed into lower dimen- 
sional operations. Consider the resampling of a digital 2D image / from a sampling pattern 
with samples Si lying on a family of curves, denoted CFi, to another another pattern with 
points S 2 lying on a different family of curves, CF 2 , producing the digital image h. If the two 
families of curves intersect at sufficient density the method that was described with reference 
to Fig. 23 for resampling from one fan to another rotated fan can be used for general curves. 
Otherwise a third family of curves CF 3 can be introduced, which intersects both CF X and 
CF 2 at a desired density. The digital image can then be resampled from CF X to its intersec- 
tions with CF 3 , then to the intersections of CF 3 with CF 2 , and finaly on CF 2 to the desired 
sampling pattern. 

This process generalizes to 3D, for example by considering surfaces instead of curves, 
to first reduce the process to one of resampling on surfaces, and then using the 2D process to 
further reduce resampling on surfaces to resampling on curves. 

Digital coordinate transformations and resampling can often be combined, improving 
the computational efficiency. For example, the digital x-shear-scale operations used in the 
shear-scale algorithm shown in Fig. 13 can be decomposed into resampling operations on 
individual horizontal lines. 



1 Divergent-beam Fast Hierarchical Backprojection Algo- 
rithms 

1.1 Fan-beam projection and backprojection 

Consider the case of equiangular-spaced detectors located on a circle around the ob- 
ject. The fan-beam tomographic projection, at a source-angle p, of a two-dimensional image 
f(x, y) is denoted by {Hpf)(i) and is defined as the set of line integrals along the rays of 
a fan, parameterized by 7, centered at the source position on a circle of radius D from the 
origin. The function f{x) is assumed to be zero-valued outside a disc of radius R. 
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The fan-beam ray V(fi, 7, T), shown in Fig. 17, is defined by 

V 05, 7, T) = D [ - ] + T [ JSJg^ ] (41) 

The fan-beam projection at source-angle /3 and fan-angle 7 is (Top/) (7) = Jrlo f( v iP> 7. r))dT. 
Since / is zero outside the disc of radius R, the integral needs only be performed within the 
disc i.e. between Tst and Tend- 

In computed tomography with the fan-beam geometry, projections are available at a 
set of discrete source angles {p p : p = 1, 2, P}, and within each fan the angles of the 
rays are indexed by {7,- : j = 1, 2, J}. In the case of equiangular fan-beam geometry 
the detectors are equally distributed on the arc of a circle centered at the source, so the 
fan-angles are evenly spaced. In the case of equispaced detectors, the detectors are equally 
distributed on a line perpendicular to the line from the source position to the origin. The 
fast backprojection algorithms for fan-beam geometry described here assume equiangular 
distribution with J odd and 7j = A 7 • (j - ( J + l)/2). However, the algorithms may be 
easily extended to other fan-angle geometries. 

The reconstruction algorithm from a set of P fan-beam projections {R Pi (7) : i = 
1, 2, P}, may be expressed as the scaling and filtering of each fan-beam projection fol- 
lowed by a weighted backprojection (5): 

/(*) = £ W(T(£, p P )) qf ,M*> AO) W 
p=l 

where W(T) is an appropriate weight function, T((a>), 0)) and nf((x),P) are the distance 
along the ray between source and image point x for source position p, and the fan angle of 
that ray, respectively, and qp(j) are the weighted and filtered projections 

Q fi (j) = #0(7) * 5(7), R'ed) = Mnr) • D • cos 7 

where g(j) is an appropriate filter. 

Definition 1.1 The weighted backprojection of a function q at a single source angle 0 is 
defined by 

(W^)(x) 4 W{T{x, 0)) 9 (7(*, P)) ( 43 ) 
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and the weighted backprojection operator W$ (where 0 € [0, 2n) p ) is defined by 

W*} W (*) = E(W^ p9 p) (S ) (44) 

P =i 

Therefore /(x) = (W^{? p }^ =1 )(x). 
It is easily shown that 

Wp = /C(-0)W o (45) 
where /C(£) denotes the rotation by P operator, and consequently that 

W*}5*)(*) = S (46) 

p=i 

Thus, as in the parallel-beam case, the weighted backprojection of P fan-beam projec- 
tions may be expressed as the sum of weighted-zero-backprojected images as seen in Figure 
4, with B 0 = W 0 . In fact, close analogs of the backprojection operator denned in Equa- 
tions (43) and (44) apply more generally to the reconstruction of functions in two and higher 
dimensions from projections of a general form, with appropriate definition of the functions 
7(£ , P) and T(x , P). The methods of the present invention extend to these other applications, 
with appropriate definition of a coordinate transformation K. 

* 

1.2 Fast hierarchical backprojection 

The Fast Hierarchical Backprojection Algorithm for the fan-beam geometry is similar 
to that for the parallel-beam case. It combines the fan-beam projections in a ternary hierar- 
chical structure, exploiting the fact that intermediate images formed by projections that are 
close to each other in projection angle P can be sampled sparsely. The block-diagram is 
shown in Figure 18. Similar algorithms can be derived for binary or other (possibly mixed) 
radix structures. 

The equations that govern the combination of the underlying continuous images in 
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Figure 18 are as follows : 

h,m = W o9 m (47) 

/l+l,m = /C(5,, 3 m-2)Il^m-2 + ^,3m-l + /C(5|,3m)ii,3m i = l,2,...,L (48) 

= + K(-*/2)h,2 + JC(-7r)J L , 3 + /C(-37r/2)/ L ,4 (49) 

In the embodiment described next, it is assumed that the source angles are uniformly 
spaced in [0, 2-ir) as follows : 

ft = _^(i - I) + Ap(j - 1) where = tt/(2 • 3 L ) 

The intermediate rotation angles are then chosen as in the parallel-beam case using 
Equation (22), with replacing A 0 . 

The constituent fans of an intermediate image in three levels of the algorithm with 
block diagram Fig. 18 are illustrated in Fig. 20. 

1.2.1 The fan-beam sampling scheme 

One embodiment of the fan-beam algorithm uses a sampling scheme of the intermedi- 
ate images in the algorithm derived by analogy to the parallel-beam case, which is described 
next. Alternatives and improvements are described later. 

A single-fanbeam-backprojection image (an image W 0 q produced by weighted back- 
projection of a single fanbeam ) has a structure amenable to sparse sampling. It follows from 
Equation (43)that (Wpq) (V(/3, 7, T)) = W(T)q(i). A single sample at a particular T = T x 
on the ray indexed by 7 in a fan oriented at /? = 0 is therefore sufficient to specify the value 
of the image along the entire ray: it is simply proportional to W{T) where T is the distance 
from the source to the point in the ray. 

An intermediate image in the algorithm is a sum of several rotatated versions of such 
single-fanbeam-backprojection images, each generated from a constituent projection. We re- 
fer to the fan at the angle that is at the center of the angular interval spanned by the constituent 
projections as the central constituent fan. This may correspond to an actual projection - e.g., 
in the case of a ternary algorithm, - or a to virtual projection - e.g., in the case of a binary 
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algorithm. For example, the fan in Fig. 20(a) is the central constituent fan for both Fig. 20(b) 
and Fig. 20(c). 

Recall that in the parallel-beam algorithm the intermediate images are sampled on 
cartesian sampling patterns aligned with the central constituent. In the rotation-based parallel- 
beam case the intermediate images are sampled along parallel vertical rays. The spacing of 
samples along each of these parallel rays is chosen to guarantee that the other constituent 
projections of the image are sufficiendy sampled (by the Nyquist sampling criterion) along 
that ray. The necessary spacing is exactly equal to the spacing of the intersection points of 
the vertical rays with a set of rotated parallel rays corresponding to the extremal constituent 
projection, that is, the one farthest away in view angle from the center projection. 

In direct analogy to the parallel-beam case, here the intermediate images are sampled 
along the rays of the central constituent fan. The sampling points along the rays of the 
central fan are the intersection points of the central fan with the extremal constituent fans. 
This results in samples that are not uniformly spaced along each ray. Consider two fans 
V{Pi1ji 0I/=i ^ d ViP'yjp Assuming the equiangular fanbeam geometry with an 

odd number J of detectors, 7; = A 7 - (j - ( J + l)/2). The points on the r th ray of the £ 
-fan (corresponding to fan-angle 7 r ) that intersect the £'-fan, are determined by the equation 
V (£, 7r , T) = V (/?', 7j> T) l/=i> wWch has the solution 

* 

Tp,p>,v r C7i) - gin _ p, + 7r _ ^ 

The function 2> 1/8 /, 7r (Y) describes how the fan V(p', -, •) varies along the r th ray of 
the fan V(/3, -, •). It carries information about the varying sampling rate along a ray of the 
fan. The steeper the slope, the sparser are the samples. The local sampling rate at any V is 
proportional to {dT/di)' 1 . A typical T(Y) is displayed by the dashed line in Fig. 22. 

Two modifications are made to Tp^>, lT (7') to get a new sampling function T PtP > ilT (7') : 

1. Intermediate images are to be sampled within the disc of radius R, with margins of 
at least one sample on each ray outside the disc (one sample with T < T S t and one sample 
with T > Tend)- Define 7^ = Tjj, >7r (Tend) : the value of 7' corresponding to Tend- 
Fans with adjacent source angles may lack intersection points with T > Tend- In order to 
rectify that, the slope of f p ^ nr (V) for T > Tend is clamped to that at T PtP >^ r (Y E nd)- 
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2. The final image is sampled on a Cartesian pattern with unit sampling intervals. 
This is the smallest interval at which intermediate images need to be sampled. The point 
7i at which a unit sampling rate along the j th ray is achieved is determined by solving the 
equation d/{di)Tp^ ylr {l) = V A 7 for 7. To maintain the unit sampling interval constraint, 
the slope of 2>,0', 7r (Y) for 7' > 7^ is set to 1/A 7 . 

Fig. 22 displays T p ^ t y T (j) and T^p^i) for a typical /3 and 7. It also displays the 
T-values of the actual sampling points on the ray which are the set {fp^ ilr (j A 7 ) : j e Z}. 
The integers j are chosen so that there are margins of at least one sample on each side of the 
image boundary. 

Clearly, the locations of the sample points prescribed by the principles outlined herein 
can be computed from the geometry of the fans, Equation (21) for the set N ltTn of the con- 
stituent projections, the selected rotation angles 5 li7n for the intermediate images (given in 
Equation (22) for the case of equispaced view angles), and the chosen form for T$£* ar (7). 

Separable rotation and up-interpolation for fan-beam sampling 

As described in the Overview, the upsampling operations, and upsampling combined 
with rotation operations can be decomposed into computationally efficient one dimensional 
resampling operations. 

1.2.2 Sampling Scheme Based On Local Fourier Structure 

Given an image f(x) : R n -» R, we want to find a sampling function t(x) : R n -* W 1 
such that f(t(x)) has a small essential bandwidth and therefore can be sampled on the set of 
points {t(m) : m G Z n }. We find this sampling function as follows. As will be explained 
shortly, knowing how f{x) is composed of its constituent projections, we find the matrix 
function v(x) = f^J [§ UlU ] ^ describes how the function / should be sampled locally 
at the point x. We then integrate this matrix function over the image domain to get the 
sampling function t(x). 

» 

In our algorithms, we know that the intermediate image / is of the form : / = 
J2 P fC(8 p )W 0 q p for some angles 5 P € \Pmin, Pma*]- Ignoring the weighting by W(T), and 
assuming that the projections q p are sampled exactly at the Nyquist rate, we know the lo- 
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cal structure of the spectral support at each point in the image /. Because of the fan- 
backprojection of each band-limited projection q p , an image that is fan-backprojected from 
a single projection has a negligible spatial bandwidth in the direction of the spoke of the fan, 
while it's bandwidth in the perpendicular direction is inversely proportionate to the distance 
from the vertex of the fan as seen in Figure 24. When a set of these fans are rotated and 
added together, the local spectral support of this resulting image is the union of the spectral 
supports of the individual fans. 

For example when the constituent fans have source angles between p m i n and fimax, as 
shown in Figure 25, the local spectral support at a point x in the image domain is a warped 
wedge with radii oriented between 0 min and 0 mox , and a radial bandwidth of ClR/r g at angle 
0. Here Q is the spatial bandwidth of the back-projected projection at the center of the image 
plane assuming the projections are sampled at the Nyquist rate. In the equiangular case Q is 
number of samples of the projection divided by the length of the arc through the origin of the 
image plane that the backprojected projection covers. Mathematically the spectral support is 

{(cos0u>0,sin0u;e) : 0 € [0 mini 

^and^ehni,^]} (50) 



where 



('max — taJO. ~ n ! ~~ 

R COS p max + x 2 

R Sin Pmin — Xi 

V min - tan R p ^ + x2 



Q e = RQ/rg 



r B = ^( Xl - Rsm9) 2 + {x 2 - R cos 0) 2 



Given this knowledge of the local two-dimensional fourier structure at a point x in 
the image, the matrix function v(£) at that point is the sampling matrix that, if the spectral 
support at that point were uniform across the whole image (such as in the parallel-beam 
case), would efficiently sample the image. In the intermediate images of interest here, this 
produces two distinct small-bandwidth and large-bandwidth sampling directions. We fix the 
first coordinate to be the small-bandwidth direction. 
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We approximate this sampling matrix function v(x ) by ( [ a 0 x s ° 2 ] [ I ° 1 l ] [ Q2 1 ]) • 
O min sin 9 min + Q max sin 9^ _ gmm cos fl min - ^ cos fl mc 

a2 ~~ fi min COS 9 min + "max COS 0 mo x "max Sm 0 mi „ + Onto Sin 0 mi 

(5 



Here 

max 
min 

r 1 

marc / 

(51) 



Here f2 min A fi* rofn and C2 m ax ^ Ve max • Angle 0 mia refers to the angle of the projection 
corresponding to the source-angle (0 min + A„ax)/2, and fl mid = Q 6mid - These parameters, 
chosen by geometrical arguments on the spectral support, ensure that the spectral support 
upon transformation by the sampling matrix is restricted to [-it, it] x [-it, it]. 

The sampling function for the whole image is found by integrating this sampling ma- 
trix function across the whole image — solving the set of differential equations : ^ = 
Vij (t(n)) for i, j = 1,2. This may be solved numerically. Even if this exact prescribed 
pattern is not used, the local Fourier support analysis will evaluate the effectiveness of any 
sampling pattern. 

The resulting sampling patterns for a few intermediate images are shown in Figs. 
26(A) and 26(B). For the fan-beam case, this local Fourier-based method produces sampling 
patterns that are similar to the ones resulting from the intersection-based methods described 
earlier.d The previously described separable method to resample from one sampling pattern 
to another can be used with these patterns also. This local Fourier sampling method may 
be applied directly to find sampling schemes for arbitrary projection geometries over lines, 
curves or planes over arbitrary dimensions. In the case of the parallel-beam geometry, it 
reduces to that discussed earlier. 

1.2.3 Alternative Sampling Schemes 

The sampling schemes in which the samples are located on a fan whose vertex is on the 
source trajectory (called a "sampling fan") has some shortcomings. The sampling points are 
chosen separately for each ray resulting in a scheme that is not necessarily optimal in a two 
dimensional sense. Though the final image in the algorithm is sampled on a uniform rect- 
angular grid with unit-spacing, the intermediate images using the above mentioned scheme 
are sampled more densely than necessary in certain regions of the image (eg. close to the 
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vertex of the fan). In order to rectify this the above described scheme might be modified to 
make sparse the sampling in such oversampled regions. Two such possibilities are illustrated 
in Figure 27. Both these possibilies incorporate the fact while in the first level the image is 
sampled efficiently on a fan, in the final level it needs to be sampled on a rectangular grid. 
These schemes attempt to incorporate the gradual transition from sampling on a fan to a grid. 

In Fig. 27 (A) samples are selected on a pseudo-fan whose vertex is located further 
from the origin than the source radius. In successive levels intermediate images are formed 
from fans from a larger range of source angles and the distance of the vertex of the pseudo- 
fan from the origin increases. In the final level the vertex is at infinity; i.e. the rays are the 
parallel lines of the rectangular grid. In Fig. 27 (B) the samples are located, not on a fan, 
but on a beam that is less divergent nearer the source position. In successive levels the beam 
becomes more parallel and less divergent. These or other sampling schemes that take a two 
dimensional (or, for 3D images, a 3D) point of view will contribute to a faster algorithm. 

1.2.4 Other Optimizations 

• In the fan-beam case with a full trajectory, source angles are distributed in [0, 2n). 
This allows us to take advantage of rotations by -tt/2, -7t and Stt/2 that involve no 
interpolation but only a rearranging of pixels. 

• The positions of points in the sampling patterns for all levels of the algorithm can 
be pre-computed and stored in a look up table, or computed on-the-fly along with the 
processing of data. In either case, computation and/or storage can be reduced by taking 
advantage of the smooth variation of sample position as a function of ray index, which 
may be observed in Figs. 23(A) and 23(A). 

w 

1.2.5 Oversampled Fast Hierarchical Backprojection 

As in the parallel-beam case, oversampling by a factor 7 < 1.0 can be used to improve 
accuracy. One way to achieve this in the fan-beam case, is to determine the denser sampling 
patterns using a fan-angle spacing A' e = • 7. 
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1.2.6 Algorithms for Short Scan and alternative Angle Distributions 

The methods of this invention are also applicable to the so-called short-scan data ac- 
quisition format, where the range of source angles is smaller than [0, 2ir). The difference is 
in the weighting used in backprojection: an extra weighting factor (e.g., Parker weighting) 
is introduced to account for the short scan. Likewise, the algorithms generalize to nonuni- 
form distributions of the source angle, or even small sets of source angles that may not allow 
reconstruction of a complete image. Furthermore, the order of operations in the algorithm 
may be modified in various ways (for example pipelining), which would be known to those 
skilled in the art. For example, the operations may be ordered such that processing may 
begin as soon as a single projection is available. 
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